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Hiis  is  the  final  report  of  the  project  to  apply  geostatistics  and  wavelet  analysis  to  part  of  the  area  chosen  at  Fort 
A.P.Hill  in  northeastern  Virginia.  This  rq)ort  embraces  several  facets  of  the  collaborative  work  starting  with  a 
summary  of  time  spent  at  the  US  Topographic  Engineering  Center  and  of  time  spent  by  TEC  personnel  at  the 
University  of  Reading.  An  analysis  of  cokriging  of  Korean  temperature  data  with  elevation  has  been  done  to  assess 
whether  the  accuracy  of  estimates  could  be  improved.  The  improvement  was  small  compared  with  the  tenqjerature- 
estimates  for  the  USA  done  previously.  A  comparison  between  geostatistics  and  wavelet  analysis  has  been  under¬ 
taken.  It  seems  that  wavelets  provide  a  more  accurate  method  for  data  reconstruction,  but  geostatistics  is  more  _ 
appropriate  for  exploring  different  resolutions  of  spatial  variation  that  have  been  identified  by  the  variogram.  A 
detailed  analysis  of  ground^cover  at  AP.Hill  has  been  done  for  three  surveys  (one  was  done  in  a  previous  rq)ort). 
The  results  show  a  relation  between  the  scales  of  variation  in  certain  groimd  cover  attributes  and  the  SPOT  data. 
Cross  variograms  suggest  that  there  are  spatial  relations  among  variables  and  the  image  data.  The  analysis  of 
elevation  data  shows  that  the  patterns  in  its  variation  correspond  closely  with  those  for  the  NIR  waveband.  Visbally 
the  relation  is  very  strong,  but  statistically  it  is  weaker.  '.1 
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EXECUTIVE  SUMMARY 


This  report  contains  most  of  the  two  earlier  interim  reports,  together  with  the  most 
recent  results.  In  addition,  there  are  three  Appendices;  a  published  paper,  an  abstract  for 
the  Geostatistics  Congress  in  2000,  and  some  new  computer  programs  that  have  been 
written  for  this  project.  The  results  of  applying  them  will  be  provided  in  a  possible  new 
program  of  work.  A  cokriging  analysis  of  Korean  temperature  data  with  elevation  has 
been  done  to  assess  whether  the  accuracy  of  estimates  of  temperature  could  be 
improved  using  the  elevation.  The  improvement  was  small  compared  with  the 
temperature  estimates  for  the  USA  done  previously.  The  relation  between  temperature 
and  elevation  was  strong  after  both  variables  had  been  detrended.  The  maps  show 
slight  differences  as  do  the  validation  results. 

A  major  aim  of  this  project  has  been  to  examine  the  relation  between  geostatistics  and 
wavelet  analysis  for  exploring  spatial  variation  in  imagery  at  different  spatial  scales  and 
data  reconstruction.  It  seems  that  overall  the  wavelet  analysis  provides  a  more  accurate 
method  for  data  reconstruction.  However,  it  is  not  straightforward  because  the 
reliability  of  the  restored  values  when  compared  with  the  original  data  varies  over  the 
region.  Kriging  performs  less  well  where  there  are  marked  changes  in  the  reflectance 
values  and  they  appear  to  be  non-stationaiy.  Kriging  analysis,  however,  seems  to  be 
more  appropriate  for  exploring  different  resolutions  of  spatial  variation  that  have  been 
identified  by  the  variogram.  The  variogram  could  be  used  to  make  the  wavelet  analysis 
for  different  resolutions  more  efficient  because  the  relevant  scales  could  be  targeted  at 
the  outset.  The  results  show  that  similar  patterns  of  variation  were  retrieved  by  both 
analyses  for  the  long  range/low  frequency  component. 

A  detailed  analysis  of  ground  cover  at  A.  P.  Hill  has  been  done  on  the  information  from 
three  surveys  (another  survey  was  described  in  a  previous  report).  The  results  show  a 
relation  between  the  scales  of  variation  in  certain  ground  cover  attributes  and  the  SPOT 
data:  in  particular  grass,  mixed  woodland,  forest  and  wetland.  The  multivariate 
variograms  of  the  quantitative  data  and  the  categorical  data  confirm  a  relation  between 
NIR  and  ground  cover  in  terms  of  the  spatial  scales  identified.  Cross  variograms 
between  the  ground  cover  types  and  also  between  each  type  and  each  waveband  suggest 
that  there  are  spatial  relations  among  variables  and  the  image  data. 

The  analysis  of  the  elevation  data  shows  that  the  patterns  in  its  variation  correspond 
closely  with  those  for  the  NIR  waveband.  These  results  confirm  our  earlier 
interpretation,  in  the  previous  project,  that  the  observed  changes  in  the  long-scale 
component  of  the  variation  coincides  with  changes  in  the  physiography.  The  raw  data 
and  detrended  data  were  analysed,  but  there  was  little  difference  between  the  results. 
Visually  the  relation  between  elevation  and  NIR  is  very  strong,  but  it  is  more  difficult  to 
show  this  statistically. 
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PART  I:  REPORT  ON  TIME  SPENT  WITH  TEC  PERSONNEL 


This  report  embraces  several  different  components  and  includes  the  material  in  the 
previous  interim  reports.  It  begins  with  a  brief  summary  of  the  work  that  was  done  by  Dr 
Oliver  at  TEC  in  February  1998  which  was  part  of  this  contract  (albeit  slightly  premature) 
and  July  1999,  and  the  work  done  while  E.  Bosch  and  E.  Shine  were  working  with  Dr 
Oliver  at  Reading  in  September  1998  and  May  1999,  respectively.  Part  H  of  the  report  is  a 
small  piece  of  work  for  Dr  Krause  on  cokriging.  Part  in  focuses  on  a  comparison  between 
wavelet  analysis  and  geostatistics.  Part  IV  on  the  vegetation  surveys  and  analysis  of  the 
digital  elevation  model  for  A.  P.  Hill.  There  are  three  appendices  to  the  report.  The  first  is 
a  copy  of  the  paper  presented  by  Dr  Oliver  at  the  Geoenv’98  conference  in  Valencia  last 
year,  and  which  has  subsequently  been  published  (Oliver  et  al,  1999).  The  second  is  the 
abstract  submitted  to  the  Geostatistics  Congress,  and  the  third  a  set  of  computer  programs 
written  by  Professor  Richard  Webster  to  compute  moving  variograms,  and  moving 
averages  and  variances.  These  will  be  tested  in  the  next  phase  of  the  work. 

The  work  for  the  majority  of  this  project  has  been  based  at  Fort  A.  P.  Hill  in  northeastern 
Virginia,  about  75  miles  from  Washington,  DC.  The  area  is  intensely  dissected  by  many 
small  waterways,  and  this  appears  to  have  contributed  to  the  pattern  of  variation  observed 
in  the  image. 

Report  of  visit  to  TEC  in  February  1998 

Much  of  the  first  day  at  TEC  was  spent  discussing  the  results  of  the  first  analyses  from 
Fort  A.  P.  Hill,  and  what  other  work  should  be  done.  In  addition  the  paper  that  has  now 
been  accepted  by  the  International  Journal  of  Remote  Sensing  was  also  discussed  and 
suggestions  for  improvement  considered  and  incorporated.  Since  Dr  Oliver  was  to  brief 
the  senior  management  team  at  TEC  including  Dr  Roper  the  contents  of  the  briefing  were 
also  ratified  at  the  outset.  A  further  paper  on  this  subject  was  considered  for  presentation 
at  the  Geoenv’98  conference  in  Valencia  (Oliver  et  al.  1999)  and  this  has  now  been 
published.  (Appendix  1). 

The  main  aim  of  this  visit  was  to  work  with  Jim  Shine  to  enable  him  to  make  fiall  use  of 
Genstat.  A  set  of  programs  was  prepared  to  cover  exploratory  data  analysis  (histograms, 
box  plots,  summary  statistics,  trend  detection  and  so  on),  variogram  analysis  and  modelling 
and  kriging.  Al  of  the  programs  were  examined.  They  were  eventually  compiled  as  part 
of  the  aide  memoire  that  formed  an  Appendix  of  the  final  report  for  the  previous  contract 
(Contract  N68171-97-C-9029). 

Jim  Shine  and  Dr  Oliver  worked  though  all  of  the  programs.  A  problem  was  identified 
with  the  kriging  algorithm  in  Genstat  which  was  eventually  reported  to  the  NAG  library 
and  corrected.  TEC  then  received  a  new  implementation  of  the  package.  At  least  half  of 
the  time  at  TEC  was  spent  instructing  Jim  in  the  use  of  the  programs  and  interpreting  the 
results.  In  addition  we  had  several  discussions  on  geostatistics. 

Haifa  day  was  spent  on  the  briefing  to  Dr  Roper  and  senior  staff  at  TEC,  and  in  answering 
questions  arising  from  this.  During  the  course  of  our  collaboration  we  have  covered  a 
substantial  amount  of  work  and  much  of  it  was  described  briefly  at  this  meeting.  Dr  Roper 
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showed  considerable  interest  in  what  has  been  done  and  when  he  visited  the  University  of 
Reading  in  November  1998  it  was  clear  that  he  had  a  sound  appreciation  of  the  value  of 
geostatistical  analysis.  The  discussion  that  followed  the  briefing  was  stimulating  and  well 
considered. 

Other  discussions  were  held  with  Edward  Bosch  about  comparing  wavelets  and 
geostatistics.  This  culminated  with  an  arrangement  for  him  to  visit  Reading  in  September 
1998. 


Visit  by  E.  Bosch  to  University  of  Reading  September  1998 

Dr  Oliver  and  E.  Bosch  worked  together  for  a  week.  The  time  was  used  for  analyses, 
interpreting  results  and  discussion.  Several  analyses  were  undertaken  -  some  of  which 
feature  in  the  report.  Others  have  been  done  by  both  of  us  subsequently.  The  visit  was  very 
profitable  to  both  of  us.  As  a  result  of  this  investigation  we  have  submitted  an  abstract  to 
the  Geostatistics  Congress  to  be  held  in  April  2000,  and  this  has  now  been  accepted  for 
presentation  and  publication.  The  appendix  is  appended  at  the  end  of  the  report.  The  paper 
will  acknowledge  the  support  of  US  Army  and  of  TEC  in  this  work,  and  will  be  authored 
jointly  by  M.  A.  Oliver,  E.  Bosch  and  K.  Slocum, 


Visit  by  James  Shine  to  University  of  Reading  May  1999 


Dr  Oliver  and  Mr  Shine  worked  together  for  a  week  in  May  1999  when  Mr  Shine  visited 
the  University  of  Reading.  This  time  was  used  for  analyses,  a  draft  outline  of  a  proposed 
paper  and  discussion.  Mr  Shine  wished  to  go  over  the  analysis  for  computing  the 
variogram  from  large  sets  of  data.  We  experimented  with  some  of  the  1-m  data  for  A.  P. 
Hill  using  the  program  ggridS.f,  written  for  the  project  by  Professor  R.  Webster.  Mr  Shine 
wanted  to  develop  his  experience  in  this  so  that  he  can  compute  variograms  from  large  data 
sets  within  a  short  time.  He  left  reading  feeling  confident  about  this.  In  addition  we  also 
fitted  models  to  the  variograms  with  Genstat  and  again  this  reinforced  what  we  did 
together  at  TEC  last  year. 

A  considerable  part  of  the  week  was  spent  discussing  the  results  from  the  final  report  of 
contract  No.  N68171-97-C-9029  which  we  now  wish  to  publish.  We  examined  previous 
issues  of  the  International  Journal  of  Remote  Sensing  to  see  whether  this  was  suitable  for 
this  work.  We  decided  that  it  was,  but  that  as  the  content  will  be  small  compared  with  the 
previous  paper  we  shall  submit  it  as  a  Letter.  This  is  confusing  because  this  form  of 
publication  is  a  short  paper  in  essence  and  will  suit  our  needs  perfectly  in  this  instance.  An 
outline  of  the  paper  has  been  prepared  and  the  introduction  written.  We  shall  continue  with 
this  when  Dr  Oliver  visits  TEC  in  July. 

The  remaining  time  was  spent  discussing  the  recent  work  on  the  ground  survey  data.  Part 
of  this  work  is  included  in  this  report.  However,  there  is  still  some  way  to  go  on  this.  We 
also  discussed  future  work.  One  idea  is  to  compute  a  moving  variogram  to  deal  with  the 
problems  of  local  trends  or  non-stationarity  in  the  data.  This  arises  at  A.  P.  Hill  for 
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example  where  there  are  water  bodies  and  areas  of  hard  standing  and  buildings.  The 
computer  code  for  this  will  be  written  as  part  of  the  current  contract,  but  any  testing  of  it 
will  have  to  be  done  in  the  future. 


Report  of  visit  to  TEC  in  July  1999 

Dr  Oliver  visited  TEC  in  July  1999  for  three  days.  On  arriving  she  gave  a  short  briefing  to 
Mr  W.  Clarke  (head  of  section)  on  the  status  of  our  current  research,  how  this  builds  on 
work  done  in  the  past  and  where  any  future  research  is  likely  to  develop.  On  the  second 
day  Dr  Oliver  had  a  meeting  with  Dr  Roper  together  with  Mr  Shine.  This  was  to  discuss 
present  work  and  also  spatial  investigations  more  generally.  Dr  Roper  invited  Dr  Oliver  to 
give  a  general  briefing  to  TEC  next  year  on  the  research  to  date. 


Part  of  each  day  was  spent  with  Mr  E.  Bosch.  We  have  been  exploring  a  one-dimensional 
set  of  radon  values  in  soil  where  we  know  there  are  distinct  boundaries.  The  aim  is  to  see 
how  wavelet  analysis  deals  with  this  variation  and  also  that  of  the  residuals  from  the 
geological  classes.  We  explored  different  levels  of  resolution  for  the  raw  data.  This  work  is 
still  to  be  completed. 

The  work  with  Mr  Shine  began  by  extracting  part  of  the  data  from  the  SPOT  image  and  the 
digital  elevation  model  (DEM).  We  plan  to  explore  the  relations  in  this  smaller  file  in  more 
detail  because  statistically  the  relation  between  the  wavebands  and  the  DEM  was  weak,  yet 
it  was  fairly  strong  for  the  NIR  band  visually.  The  weak  relation  might  arise  from  the  areas 
of  hard  standing  and  buildings  which  have  no  particular  relation  with  the  elevation.  The 
program  ggrM/ would  not  work  with  these  small  files  -  Mr  Shine  has  since  discovered  that 
the  zero  origin  has  caused  part  of  the  problem. 

We  continued  the  discussion  about  the  Letter  for  IJRS  and  have  decided  to  use  NDVI  of 
subsets  from  the  whole  site  covered  by  the  1  m  data.  The  additional  work  to  prepare  the 
new  variograms  for  this  has  now  been  completed. 
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Part  n 

Cokriging  temperature  data  in  Korea 

The  data  for  the  analysis  were  provided  by  Dr  P.  Krause.  They  comprised  temperature  and 
elevation  records  at  100  sites  irregularly  scattered  over  Korea.  In  addition  elevation  had 
been  measured  at  another  565  sites.  Table  1  gives  the  summary  statistics  for  these  variables 
at  places  where  they  were  both  measured.  Both  have  distributions  that  depart  from 
normality,  in  particular.  Although  a  geostatistical  analysis  does  not  assume  that  the  data  are 
normally  distributed  it  is  generally  advisable  to  transform  the  data  to  a  near-normal 
distribution  for  the  variogram  analysis  to  stabilize  the  variances. 

Both  variables  were  transformed  to  common  logarithms  and  for  elevation  the  skewness 
decreased  markedly  and  the  transformed  data  are  close  to  normal.  Temperature  departs  less 
so  from  a  normal  distribution,  but  after  transformation  to  common  logarithms  the  departure 
from  normality  increases. 


Table  1:  Summary  statistics  for  Elevation  and  Temperature 


Elevation 

Temperature 

Log  Elevation 

Log  Temperature 

Number  of  observations 

100 

100 

100 

100 

Mean 

403.45 

53.02 

5.17 

3.97 

Minimum 

8.00 

33.00 

2.08 

3.50 

Maximum 

4546.00 

62.00 

8.42 

4.13 

Variance 

574928.23 

24.95 

1.53 

0.011 

Standard  deviation 

758.24 

4.99 

1.24 

0.103 

Skewness 

3.836 

-1.45 

0.21 

-1.93 

The  data  were  also  examined  for  trend  as  part  of  the  exploratory  data  analysis.  This  would 
generally  be  normal  practice  when  one  of  the  variables  is  elevation  because  it  can  vary  in  a 
predictable  way.  However,  in  this  case  it  was  temperature  not  elevation  whose  variation 
comprised  a  large  element  of  trend.  For  elevation  linear  trend  counted  for  13.8%  of  the 
variation,  and  quadratic  trend  for  21.0%.  This  is  much  less  than  expected.  It  is  marginal  as 
to  whether  this  degree  of  trend  should  be  removed,  but  it  was  to  ensure  that  the  analysis 
was  reliable.  For  temperature  the  trend  was  much  greater:  a  linear  trend  accounted  for 
74.9%  of  the  variation  and  the  quadratic  one  77.9%.  Clearly  a  linear  trend  model  is 
adequate  for  describing  the  trend  for  temperature. 

The  aim  of  this  analysis  was  to  assess  whether  temperature  could  be  estimated  more  reliably 
with  the  use  of  additional  information  from  elevation.  In  geostatistics  the  method  used  is 
known  as  cokriging.  The  value  of  the  method  is  that  it  can  be  used  to  estimate  a  property 
that  is  more  expensive  to  measure  using  information  from  another  variable  with  which  it  is 
coregionalized  and  that  is  cheaper  to  measure  or  that  does  not  change  with  time.  This  is 
particularly  true  in  general  for  temperature  and  elevation.  There  is  a  physical  reason  for  their 
relation  and  elevation  does  not  change  substantially  in  the  short  term.  Therefore,  once  a 
digital  elevation  model  has  been  produced  it  is  a  source  of  inexpensive  and  reliable 
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information.  Cokriging  depends  on  the  two  (or  more)  variables  being  strongly  coirelated. 
From  the  correlation  matrix  below  it  is  clear  that  the  correlation  between  elevation  and 
temperature  is  moderate. 


Table  2:  Correlation  matrix  for  temperature  and  elevation  in  Korea. 

***  Correlation  matrix  *** 

Elevation  1  1.000 

Temperature  2  -0.741  1.000 

1  2 

This  level  of  correlation  would  suggest  that  it  is  worthwhile  pursuing  a  coregionalization 
analysis.  The  classical  correlation  coefficient  does  not  take  spatial  location  into  account, 
therefore  the  relation  spatially  could  be  either  better  or  worse. 


Cokriging:  Theory 

The  cross  variogram 

This  is  the  logical  extension  of  ordinary  kriging  to  situations  where  two  or  more  variables 
are  spatially  interdependent  or  co-regionalized.  The  first  stage  is  to  model  the 
coregionalization.  The  two  regionalized  variables,  Za(x)  and  Zv(x),  denoted  by  u  and  v, 
both  have  an  autovariogram  defined  by: 

/„(h)  =  |E[{Z„(x)-Z„(x  +  h)}^] 

and 

j.,(h)  =  |E[{Z,(x)-Z„(x  +  h)}^], 

and  a  cross  variogram  defined  as: 

r^(h)  =  |e[{Z„  (x)  -  Z„  (x  +  h)}  {Z,  (x)  -  Z,(x  +  h)}]. 

The  cross  variogram  function  describes  the  way  in  which  w  is  related  spatially  to  v. 
Provided  that  there  are  sites  where  both  properties  have  been  measured  yav(h)  can  be 
estimated  by: 

1  m(h) 

(h)  =  Z  [{Za  (x)  -  Za  (x  +  h)}  {2v  (*)  "  (*  +  h))  ]  • 

which  provides  the  experimental  cross  variogram  for  u  and  v. 
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The  cross  variogram  can  be  modelled  in  the  same  way  as  the  autovariogram,  based  on  the 
linear  model  of  coregionalization.  Each  variable  is  assumed  to  be  a  linear  sum  of 
orthogonal  random  variables  Y (x); 

t=i  j=i 

in  which 
E[Z„(x)]  =  /4,, 

(x)  -  T,*  (X  +  h)}  {F/  (x)  -  Y^'  (x  +  h)}] 

=  gi^  (h),  positive  for  k  -  k'  mdj  =  f 
=  0  otherwise 

The  variogram  for  any  pair  of  variables  u  and  v  is; 

Jfc=l  J=l 

We  can  replace  the  products  in  the  second  summation  by  to  obtain: 


The  i*^are  the  nugget  and  sill  variances  of  the  independent  components  if  they  are 
bounded,  and  for  unbounded  models  they  are  the  nugget  variances  and  gradients. 


Cokriging 

Once  to  coregionalization  has  been  modelled  it  can  be  used  to  predict  the  spatial  relations 
between  two  or  more  variables  by  cokriging.  There  are  generally  two  reasons  for  using 
cokriging: 

1.  Where  one  variable  is  under-sampled  compared  with  another  with  which  it  is 
correlated.  The  sparsely  sampled  property  can  be  estimated  with  greater 
precision  by  co-kriging  because  the  spatial  information  from  the  more  intensely 
measured  one  is  used  in  the  estimation.  The  increase  in  precision  depends  on  the 
degree  of  under-sampling  and  the  strength  of  the  coregionalization. 

2.  When  values  of  all  of  the  variables  are  known  at  all  sample  points,  cokriging  can 
improve  the  coherence  between  the  estimated  values  by  taking  account  of  the 
relation  between  them. 
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If  there  are  V  variables,  /  =  1,2,...,  V,  and  the  one  to  be  predicted  is  u,  which  in  our  case 
has  been  less  densely  sampled  than  the  others.  In  ordinary  cokriging  the  estimate  is  the 
linear  sum: 

/=i  j=i 

where  the  subscript  I  refers  to  the  variables,  of  which  there  are  V,  and  the  subscript  /  refers 
to  the  sites,  of  which  there  are  «/  where  the  variable  /  has  been  measured.  The 
Xu  are  the  weights,  satisfying: 

ni  ni 

2  =  1,  /  =  M ;  and  I  ^u. 

(=1  i=l 

These  are  the  non-bias  conditions,  and  subject  to  them  the  estimation  variance  of 
Z^(B)  for  a  block,  B,  is  minimized  by  solving  equations  : 

r  ni 

for  all  v=l,2  to  V and  ally=l,2  to 

l=l  i=l 

The  quantity  yiv(jc,,  Xj)  is  the  cross  semivariance  between  variables  /  and  v  at  sites  i  and  j, 
separated  by  the  vector  f^(Xj,B)  is  the  average  cross  semivariance  between  a  site  j 

and  the  block  B,  and  \|/v  is  the  Lagrange  multiplier  for  the  vth  variable.  The  cokriging 
variance  is  obtained  from: 


V  n, 


;=i  i=i 


where  (B,B)  is  the  integral  of  (h)  over  B,  i.e.  the  within-block  variance  of  u. 


Analysis  and  results  of  cokriging 

Cross  variogram 

The  experimental  autovariograms  for  the  raw  values  of  elevation  and  temperature  were 
computed  first.  They  showed  some  similarity  in  their  shapes  and  also  ranges  of  spatial 
dependence  (Figure  1).  The  autovariograms  were  then  computed  on  the  residuals  from  the 
linear  trend  for  temperature  and  on  the  residuals  from  the  quadratic  trend  for  elevation.  In 
addition  the  elevation  was  transformed  to  common  logarithms  and  the  variogram  was  also 
computed  from  the  transformed  data.  Considering  that  the  level  of  skewness  is  substantial 
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reducing  it  appears  to  have  had  little  effect  on  the  variogram.  In  fact  it  is  less  clearly 
bounded  and  less  related  to  the  variogram  of  temperature  than  that  for  the  raw  data.  The 
variograms  computed  from  the  residuals  were  more  erratic  and  more  difficult  to  model 
than  those  of  the  raw  data.  Since  the  trend  appears  to  be  regional  in  the  case  of 
temperature,  at  the  longer  lags,  I  decided  to  do  the  analysis  on  the  raw  data  and  the 
residuals.  For  kriging  it  is  the  first  few  lags  that  are  important  and  these  are  less  likely  to  be 
affected  by  the  trend  than  the  longer  lags. 

Although  it  is  important  to  check  the  data  in  this  way,  the  changes  did  not  appear  to 
improve  the  variogram  substantially.  This  will  become  evident  when  the  cokriging  results 
are  discussed.  However,  cokriging  was  carried  out  on  the  raw  data  and  the  detrended  data. 
During  the  remaining  time  on  the  project  I  might  do  some  further  tests,  but  I  do  not  expect 
any  major  changes. 

The  experimental  auto-  and  cross-variograms  for  the  raw  data  are  given  in  Figurel.  They 
have  a  similar  form  and  the  individual  autovariograms  were  fitted  best  by  an  exponential 
model  with  a  distance  parameter  of  about  0.86  units.  The  same  form  of  model  must  fit  all 
of  the  variograms  and  the  range  or  distance  parameter  must  be  the  same.  The  nugget 
variance,  the  sills  of  bounded  models  and  the  slope  of  unbounded  models  can  be  different. 
The  coregionalization  was  modelled  by  an  exponential  function  with  a  distance  parameter 
of  0.86  units  of  latitude  and  the  lower  triangle  of  the  sills  is  given  below.  The 
coregionalization  of  the  residuals  for  elevation  and  temperature  were  also  modelled  and  the 
values  used  for  kriging.  The  variograms  for  the  residuals  were  fitted  best  by  a  spherical 
function  with  a  range  of  1 . 0 1  units  of  latitude. 


Table  3:  Models  of  coregionalization  fitted  to  the  raw  data  and  the  residuals  from  the 
trend  for  temperature  and  elevation. 

Fitted  sills  in  lower  triangle  for  the  raw  data  Fitted  sills  in  lower  triangle  for 

the  residuals 


Nugget 

Variances 

0.0 

0.0 

0.6 

Elevation 

Cross  Temperature 

0.0 

0.0 

1.4 

Nugget 

variances 

Sill 

Variances 

350826.2 

-1169.6 

6.3 

Elevation 

Cross  Temperature 

258385.0 

-821.4 

2.6 

Sill 

variances 
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Temperature 


0  1  2 


Lag  distance 


Elevation 


0  1  2 
Lag  distance 


Cross  varlogram  of  temperature  and  elevation 


0  1  2 


Lag  distance 


Figure  1:  Experimental  autovariograms  of  a)  temperature  and  b)  elevation,  and  c)  the 
experimental  cross  variogram . 
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Figure  2  shows  the  experimental  cross  variograms,  the  fitted  models  together  with  the  hull 
of  perfect  correlation  (the  two  outer  lines).  The  cross  variogram  of  the  residuals  coincide 
with  the  hull  showing  a  strong  correlation.  That  for  the  raw  data  is  close  to  the  hull. 


Temp  V  Etevoibn 


b) 


Temp  ¥  [lev  residue  Is 


Figure  2:  a)  Cross  variogram  of  the  raw  data  and  b)  cross  variogram  of  the  residuals,  with 
the  hulls  of  perfect  correlation. 
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Cokriging 

The  first  analysis  was  to  test  the  modelling  and  to  assess  the  effects  on  the  estimates  of 
using  either  the  raw  data  or  the  residuals.  Twenty  five  of  the  100  sites  were  removed  from 
the  raw  data  and  the  residuals.  Using  the  models  of  coregionalization  given  above  the 
values  at  the  25  validation  points  were  estimated  by  punctual  cokriging  for  the  raw  and 
residual  data,  respectively.  In  addition  the  raw  data  were  used  for  autokriging  the 
validation  points.  The  original  values,  the  estimates  and  the  standard  errors  are  given  in 
Table  2. 

For  every  validation  point  the  cokriged  estimate  has  a  smaller  standard  error  than  the 
autokriged  estimate.  The  differences  are  small,  but  they  show  consistently  that  cokriging 
confers  a  small  benefit  in  terms  of  estimating  temperature  more  reliably.  In  addition  the 
estimates  are  consistently  closer  to  the  original  values  for  cokriging  of  the  raw  data.  For 
the  residuals  the  standard  errors  from  colaiging  are  smaller  for  15  of  the  25  validation 
points.  This  was  somewhat  surprising  in  relation  to  the  fact  that  the  variograms  of  the 
residuals  did  not  appear  to  be  an  improvement  over  that  of  the  raw  data.  For  the  residuals 
the  trend  was  added  back  so  that  the  values  could  be  compared  with  the  raw  data.  The 
estimates  are  not  as  consistently  good  as  they  are  for  cokriging  with  the  raw  data. 


Table  4:  Comparison  between  the  raw  temperature  data,  the  autokriged  estimates  and  the 
cokriged  estimates,  and  the  cokriged  estimates  for  the  residuals  and  with  the  trend  added 
back. 

Original  Autokriging  Cokriging  Cokriging  residuals 


X 

Y  Value  Estiiuate  SE  Estimate 

SE  Estimate  Est+trend  SE 

-127.05 

37.90 

54.0 

53.33 

2.48 

53.32 

2.43  0.5392 

53.48 

2.40 

-127.10 

37.70 

54.0 

54.35 

1.73 

54.25 

1.67  0.8142 

53.71 

1.91 

-126.50 

33.50 

60.0 

60.03 

1.63 

60.07 

1.56  -1.3596 

58,38 

1.82 

-128.10 

35.20 

57.0 

57.19 

2.19 

57.36 

2.12  -0.0677 

57.96 

2.16 

-127.75 

37.90 

53.0 

53.86 

1.55 

53.64 

1.49  0.5506 

53.74 

1.79 

-128.00 

36.20 

54.0 

55.17 

4.09 

55.13 

4.08  -0.4151 

56.08 

3.84 

-126.60 

37.50 

54.0 

53.17 

2.63 

53.24 

2.60  0.3207 

54.11 

2.57 

-128.90 

37.10 

48.0 

54.50 

4.52 

54.47 

4.51  -0.0800 

55.76 

4.03 

-129.40 

37.00 

55.0 

55.01 

5.43 

55.01 

5.42  0.1675 

56.91 

4.31 

-126.75 

34.30 

58.0 

58.36 

5.29 

58.33 

5.28  -0.5767 

58.26 

4.30 

-127.65 

37.45 

56.0 

53.82 

2.99 

53.75 

2.97  0.2831 

54.34 

2.84 

-125.65 

39.60 

50.0 

49.56 

4.68 

49.63 

4.66  0.2680 

49.52 

4.02 

-129.01 

35.10 

59.0 

57.95 

1.96 

57.87 

1.93  -0.5659 

58.47 

2.09 

-124.80 

40.45 

49.0 

49.22 

4.73 

49.36 

4.72  0.6924 

48.36 

3.86 

-128.30 

41.80 

33.0 

42.68 

4.71 

42.71 

4.68-2,4696 

40.90 

3,72 

-128.60 

35.90 

57.0 

56.81 

1.71 

56.80 

1.63  -0.0015 

57.47 

1.83 

-126.50 

36.75 

54.0 

53.60 

3.01 

53.22 

2.97  -1.7979 

53.47 

2.74 

-127.10 

37.45 

54.0 

54,94 

1.93 

54.89 

1.89  1.0006 

54.88 

2.10 

-128.20 

36.40 

58.0 

54.79 

3.35 

54,76 

3.33  -0,4044 

55.91 

3.15 

-127.95 

37.40 

53.0 

53.25 

0.98 

53.30 

0.93  0.1380 

54.47 

1.53 

-129.40 

36.03 

58.0 

56.87 

1.50 

56.92 

1.43  0.5288 

58.85 

1.69 

-124.65 

38.00 

52.0 

51.97 

1.66 

51.68 

1.60  -0.4415 

53.86 

1.81 

-126.40 

34.80 

58.0 

57.45 

5.02 

57.54 

5.00  -0.6245 

57.73 

3.94 

-125.80 

39.25 

51.0 

49.95 

4.83 

50.00 

4.82  0.1461 

50.22 

4.18 

-130.40 

42.30 

45.0 

47.26 

6.34 

47.38 

6.33  -0.3178 

45,31 

4.23 
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The  entire  data  set  was  cokriged  as  above,  but  this  time  using  all  of  the  elevation  data.  The 
estimates  and  the  standard  errors  were  mapped.  Figures  3  to  5.  Figures  3a  and  4a  show  the 
maps  of  temperature  from  autokriged  and  cokriged  estimates,  respectively.  There  is 
remarkably  little  difference  between  them.  Figure  5a  shows  the  results  of  cokriging  using 
the  residuals  and  then  adding  the  trend  back.  This  is  more  different.  This  appears  to  show 
some  distortion,  however,  it  is  difficult  to  be  certain  because  we  did  not  have  the  outline  of 
Korea  to  superimpose  on  the  estimates.  This  will  be  done  at  TEC.  Figures  3b,  4b  and  6b 
show  the  standard  errors  for  temperature.  They  are  slightly  less  for  cokriging.  These  values 
show  the  pattern  of  sampling  and  also  the  coastline  of  the  country. 


Figure  3:  a)  Map  of  estimates  from  autokriging  of  temperature  for  Korea, 
b)  map  of  the  standard  errors  from  autokriging  of  temperature 
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Cokriged  estimates  of  residuals  of  temperature  with  trend  added  back  for  Korea 


62.00 

60.00 

58.00 

56.00 

54.00 

52.00 

50.00 

48.00 

46.00 

44.00 

42.00 

40.00 

38.00 

36.00 


Standard  errors  of  residuals  from  cokriging  temperature  for  Korea 


-131.00  -130.00  -129.00  -128.00  -127.00  -126.00  -125.00  -124.00 


Figure  5:  a)  Map  of  the  cokriged  estimates  of  the  residuals  for  temperature  with  trend  added  back  for  Korea, 
b)  map  of  the  standard  errors  from  cokriging  the  residuals  of  temperature. 
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PART  ni:  Comparing  wavelets  and  kriging  for  exploring  nested  scales  of  variation 

M.  A.  Oliver  and  E.  Bosch 


This  work  has  continued  to  use  part  of  the  SPOT  image  around  Anderson  Camp.  However, 
it  is  slightly  different  in  extent  from  that  used  previously  to  accommodate  the  wavelet 
analysis.  To  apply  the  discrete  wavelet  transform  we  need  the  rows  and  columns  of  the 
area  to  be  analysed  to  be  some  power  of  2;  the  area  chosen  was  2’  x  2’  =  128  x  128.  As  a 
result  the  variogram  analysis  and  factorial  kriging  had  to  be  redone  so  that  the  results  relate 
to  this  new  region  and  are  comparable  with  the  subsequent  wavelet  analyses.  The  theory 
of  factorial  kriging  was  given  in  the  final  report  for  Project  3  (Contract  number  N68171- 
97-C-9029).  As  an  addendum  to  this  report  the  paper  presented  at  the  geoENV98 
conference  is  appended  and  this  includes  the  theory  of  factorial  kriging.  The  application 
and  theory  of  wavelets  is  only  summarised  briefly  here. 

Summary  of  wavelet  analysis 

Wavelets  have  some  similarities  with  the  windowed  Fourier  functions.  Both  of  them  have 
their  energy  well  localized  in  time.  This  means  that  these  functions  decay  rapidly  in  time 
or  space,  i.e.  they  go  to  zero  fast,  throughout  the  whole  real/complex  line.  However,  they 
can  also  have  compact  support,  they  are  non-zero  in  a  finite  interval.  In  addition,  the 
respective  Fourier  transforms  of  these  functions  have  their  energy  concentrated  about  a 
small  set  of  frequencies.  An  advantage  of  the  wavelet  transform  over  the  windowed 
Fourier  transform  is  that  low  frequency  and  high  frequency  resolutions  can  be 
characterized  simultaneously.  This  means  that  wavelet  analysis  is  suitable  for  situations 
where  there  are  different  levels  of  resolution  of  variation  superimposed  on  each  other 
(Daubechies,  1992).  Wavelets  are  also  good  for  describing  transient  data,  whereas  the 
Fourier  transform  is  not.  Wavelet  analysis  is  not  affected  by  local  non-stationarity  and  this 
is  an  advantage  it  has  compared  with  geostatistics,  which  assumes  that  the  data  are  at  least 
quasi-stationary  (i.e.  locally  stationary).  Local  non-stationarity  can  arise  where  there  are 
marked  boundaries  that  result  in  a  marked  change  in  the  local  means  of  the  variable  of 
interest  (see  Part  V  of  this  report). 

Wavelets  are  oscillatory  components  that  operate  locally.  The  wavelet  analysis  starts  with 
the  choice  of  a  mother  wavelet,  w(0,  which  is  fixed.  The  mother  wavelet  can  be  dilated  or 
shrunk  to  examine  components  in  the  variation  that  occur  at  different  spatial  or  temporal 
scales.  This  enables  multi-resolution  analysis  where  different  levels  of  variation  are 
superimposed  on  one  another  (Mallat,  1998).  This  is  our  first  aim  in  this  investigation. 
The  second  relates  to  redundancy,  which  is  a  major  problem  with  image  data  because  of 
the  amount  of  information  involved.  Wavelets  are  also  of  great  value  for  data  compression 
because  they  are  able  to  remove  redundant  information  and  to  retain  the  important 
structure  of  the  data. 


Theory 

Wavelet  analysis  allows  a  signal  (information)  to  be  represented  in  terms  of  a  set  of  basis 
functions,  i.e.  basis  vectors  or  kernels.  The  basis  functions  are  a  set  of  linearly 
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independent  functions  that  can  be  used  to  produce  all  admissible  functions  of  J{t)  (Strang 
and  Nguyen,  1996).  Choosing  the  basis  functions  determines  the  kind  of  information  that 
can  be  extracted.  Thus,  a  function^/)  can  be  expressed  as 

f(t)=^^bsuWsi:(t)  (1) 


where  bsu  is  the  wavelet  coefficient  and  Wsu(t)  is  the  wavelet  at  scale  s,  translated  by  u.  A 
special  feature  of  the  wavelet  basis  is  that  its  elements  H’ja(t)  are  constructed  by  scaling  and 
translating  a  single  mother  wavelet  w(/): 

w^{t)=w(tlT -u)IT'^  (2) 

Using  Fourier  techniques,  Mallat  (1998)  shows  how  to  construct  a  wavelet  function  w(t). 
The  wavelet  function  w(t)  is  constructed  from  a  scaling  function  ^t)  which  can  be 
specified  by  a  discrete  filter  l[n].  This  means  that  you  can  obtain  the  discrete  filter  l[n] 
from  the  scaling  fimction  or  you  can  obtain  the  continuous  function  jS(t)  from  the 
discrete  filter  l[n].  The  proofs  are  not  trivial  but  he  provides  a  thorough  explanation. 
Associated  with  the  wavelet  function  w(t)  is  a  discrete  filter  h[n].  The  filter  l[n],  associated 
with  the  scaling  function  ^t),  is  a  low  pass  frequency  filter  while  the  filter  h[n] ,  associated 
with  the  wavelet  w(t),  is  a  high  pass  frequency  filter.  The  low  pass  filter  smoothes  the 
signal  while  the  high  pass  filter  retains  the  detail  in  the  signal.  This  is  no  coincidence  since 
these  are  the  properties  that  are  sought  in  the  construction  of  such  filters  and  functions. 
Furthermore,  it  is  these  filters,  and  not  the  continuous  functions  themselves,  which  are  used 
to  compute  the  Discrete  Wavelet  Transform.  That  is,  the  results  obtained  by  convolving  a 
signal  with  these  filters  are  the  same  as  those  obtained  by  convolving  the  signal  with  the 
continuous  functions.  This  is  a  remarkable  feature  since  the  filters  consist  of  a  smaller 
number  of  elements,  thus  reducing  the  amount  of  computation. 

As  was  mentioned  above,  some  wavelets  have  compact  support.  Mallat  (1998)  shows  that 
the  scaling  function  ^t)  has  compact  support  if  and  only  if  the  discrete  filter  l[n]  has 
compact  support.  Furthermore,  their  support  is  equal.  Also,  if  the  support  of  and  l[n] 
is  [Ni,N2],  then  the  support  of  the  wavelet  function  w(t)  is  [(Ni-N2+l)/2,(N2-Ni+iy2]. 

Wavelets  can  be  orthogonal  which  implies  that 

This  means  that  the  above  integral  will  be  zero  when  s  ^  S  or  u  ^  U.  That  is,  (3)  will  be 
zero  at  different  scales  5  or  translates  u  of  w(t).  Also,  when  s-S  and  u  -  U,  then  (3)  is 
equal  to  1.  This  leads  to  a  simple  formula  for  each  coefficient  isi/in  the  expansion  ofy(0- 
The  expansion  in  equation  (1)  is  multiplied  by  wjk  and  integrated  by; 


=  b,^^Jw^(t)fdt.  (4) 

All  other  terms  disappear  because  of  the  orthogonality.  Mallat  (1998)  shows  how  the 
components  bsu  are  computed  with  the  discrete  filters  l[n]  and  h[n] . 
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A  wavelet  needs  to  satisfy  the  following  condition 

+00 

\w{t)dt  =  Q  (5) 

—00 

for  perfect  reconstruction.  That  is,  equation  (5)  needs  to  be  satisfied  to  be  able  to 
reconstruct  exactly  f(t)  from  the  Forward  Wavelet  Transform. 


Multiresolution 

The  dilation  of  the  mother  wavelet  allows  us  to  analyse  the  signal  at  different  levels  of 

resolution.  Again,  instead  of  using  the  scaling  function  and  the  wavelet  function  to 
compute  the  discrete  wavelet  transform  of  a  discrete  signal  _/(«),  we  use  the  two  filters 
described  above,  low  pass  l[n]  and  high  pass  h[n] .  We  assume  that  the  number  of  samples 
of  the  signal  f(n)  is  a  power  of  two.  To  obtain  the  first  level  of  the  multiresolution 
decomposition  of  the  signal we  do  as  follows: 

We  apply  the  high  frequency  filter  h[n]  to  the  signal  f(n).  This  portion  of  the  transform 
contains  the  fine  detail  structure  of  the  signal.  That  is,  h[n]  filters  out  the  smooth 
segments  of  the  signal  and  retains  the  sharp  transitions  or  discontinuities  of  the  signal.  The 
dilations  and  translations  of  the  continuous  wavelet  transform,  are  reflected  in  the  discrete 
filters  l[n]  and  h[n]  as  decimation  and  translation,  respectively.  Thus,  each  filter  produces 
half  the  number  of  sample  points  as  those  of  f(n).  Since  we  are  dealing  with  a  finite 
number  of  sampling  points,  we  can  express  the  application  of  the  filter  h[n]  to  the  signal /a 
=f(n)  in  terms  of  matrix  and  vector  notation  fhi=  'R\fo,  where  Hi  is  a  matrix  whose  rows 
comprise  the  elements  of  the  filter  h[n]  and  fo  is  the  vector  containing  the  elements  o^f(n). 
When  this  has  been  done,  we  then  apply  the  low  frequency  filter  l[n]  to  the  signal  f(n), 
whose  purpose  is  to  smooth  the  signal.  When  the  filter  is  applied  to  the  signal,  the  sample 
size  is  reduced  to  half  the  number  of  sampling  points  in  f(n).  Similarly,  we  can  express  the 
application  of  the  filter  l[n]  to  the  signal  f(n)  in  terms  of  matrix  and  vector  notation  fli  = 
Li/^a,  where  Li  is  a  matrix  whose  rows  comprise  the  elements  of  the  filter  l[n]  and  fo  is  the 
vector  containing  the  elements  of  f(n).  Note  that  with  this  notation,  H  refers  to  high 
frequency  and  L  refers  to  low  frequency.  These  two  sets  of  high  fhi  and  low/?/  frequency 
components  comprise  the  first  level  of  the  multiresolution  decomposition  of  the  signal  f(n). 
Since  these  filters  (wavelets)  are  orthogonal,  we  can  obtain  f(n)  from  fht  and  from  fli  by 
the  following: 


fo  =/=  ^t(fhi)  +  Um  =  Hl’CHi/)  +  U(hf),  (6) 

where  Hi’  and  Li’  are  the  transpose  matrices  of  Hi  and  Li  respectively.  This  means  that 

Hi’(H,/)  +  L,’(Lx/)  =  (Hi’Hi)/-+  {UUY=  (Hi’Hi  +  UUY 

=  1/=/, 


(7) 
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where  lo  is  the  identity  matrix  of  order  corresponding  to  the  size  of  the  original  signal /('wj. 

To  obtain  the  second  level  of  the  multiresolution  decomposition  of  the  signal /(^nj  we  do  as 
follows: 

The  same  two  filters  h[n]  and  l[n]  are  now  applied  to  the  reduced  resolution  signal/?/.  We 
leave  the  first  high  frequency  resolution  level  components/?!/  untouched.  Recall  that  the 
sample  size  of/?/  is  half  that  of fo.  In  matrix  notation,  when  applying  the  filter  h[n]  to/?/, 
we  have  the  high  frequency  component  of  /?/  in  fh2  =  H^/?/  and  the  low  frequency 
component  of/?/  in  fl2  =  Ljjfli.  Again,  since  the  filters  are  orthogonal,  we  can  obtain/?/ 
from  fh2  and  from fl2  by  the  following: 

flj  =  m(fh2)  +  Um  =  H2’(H2/?r)  +  UXUfli),  (8) 

which  means  that 

Ii  =  H2’H2  +  L2’L2. 

To  compute  the  subsequent  multiresolution  levels,  we  proceed  in  this  fashion  by  applying 
the  low  and  high  frequency  filters  to  the  previous  set  of  low  frequency  components.  Since 
f(n)  has  a  finite  sample  size,  we  can  only  perform  this  decomposition  a  finite  number  of 
times.  Suppose  the  signal  f(n)  has  2'*  =  16  sampling  points.  Then  this  decomposition  has  4 
steps.  The  resolution  levels  can  be  labelled  as  follows: 


[(/?,),  (fh4\  ifhsl  {fh2),  (10) 

where  ifhi)  has  8  sample  points,  (flt2)  has  4,  (fhi)  has  2  and  both  (flt4)  and  (fl4)  have  1 
sampling  point  each.  Each  ^j+jflj  contains  high  frequency  information  of  f(n)  but  at  a 
reduced  resolution  (grosser  scale).  This  analysis  allows  us  to  examine  features  locally  and 
at  different  scales.  Note  that  we  can  reconstruct  perfectly /(>ij  from 
[(flA  (fh4),  (fhsX  (fh2),  (fhi)]  since  h[n]  and  l[n]  are  orthogonal  filters. 

At  a  given  resolution,  the  scaling  functions  -  u)  form  a  basis  for  the  set  of  admissible 
signals.  The  level  is  set  by  s,  and  the  steps  at  the  level  are  1\  The  detail  at  the  level  5  is 
represented  by  the  wavelets  w(t/2'  -  u).  Multiresolution  divides  the  frequencies  into  octave 
bands,  from  w  to  2w,  rather  then  different  frequencies. 

From  this  analysis,  to  invert  from  a  grosser  resolution  to  a  better  resolution  Mallat  (1998) 
represents  this  as  follows: 

signal  at  level  s  +  1  (local  differences)  N 

+  signal  at  level  s 

details  at  level  s  +  1  (local  averages)  ^ 

Note  that  as  5  increases,  the  details  in  the  signal  are  further  degraded  or  washed  out.  Using 
this  notation  and  the  transform  coefficients  of  the  example  in  (10),  to  obtain /?3  from/?4 
and  fh^  we  do  the  following: 

(H4’/?!4)  +  (L4’/?4)  =  (H4’H4/?3)  +  {U'Lrfls)  =/?i. 
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Thus  the  signal,  such  as  the  MR  information  from  an  image,  can  analysed  at  different 
resolutions  by  the  wavelet  transform.  The  coefficients  provide  us  with  a  measure  of  the 
energy  the  basis  vector  has  at  time  t  and/or  scale  s.  The  discrete  wavelet  transform 
decomposes  the  signal  into  a  set  of  high  and  low  frequency  components  which  correspond 
to  the  coefficients  of  the  dilated  and  translated  basis  functions  ^(t)  and  ^t)  respectively. 
The  high  frequency  component  wavelet-coefficient  wj[s,u),  provides  a  measure  of  the 
variation  of  /  in  a  neighbourhood  u  whose  size  is  proportional  to  s.  This  measure  of 
variation  in  smooth  signals  is  negligible  while  the  variation  tends  to  be  significant  in 
signals  with  jump  discontinuities.  Furthermore,  increasing  the  dilation  parameter  s  while 
filtering  the  signal  produces  a  larger  region  of  integration,  which  in  turn  smoothes  the 
signal  further.  That  is,  less  detail  and  lower  frequencies  are  obtained  as  the  dilation 
parameter  increases.  On  the  other  hand,  decreasing  the  dilation  parameter  s  in  the 
convolution  process,  generates  smaller  windows  of  integration  allowing  more  detail  and 
higher  frequencies  to  come  through. 

There  are  many  different  kinds  of  wavelets  (adaptive,  continuous,  discrete,  orthogonal, 
biorthogonal,  real  and  complex),  the  most  simple  being  the  Haar.  Nevertheless,  those  of 
Daubechies  Paubechies,  1988)  have  been  used  quite  often  since  they  first  came  out. 


Analysis  of  the  A.P.  Hill  data 

SPOT  Image 
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The  part  of  the  scene  covering  Fort  A.  P.  Hill,  Figure  6,  is  slightly  smaller  than  that  used 
before  (see  report  N68171-97-C-9029),  but  it  covers  the  same  part  of  the  image.  Analyses 
were  carried  out  on  the  complete  data  set  and  on  sub-samples  of  1  pixel  in  2  for  each 
column  and  row  (or  1  pixel  from  a  block  of  4),  1  pixel  in  4  for  each  row  and  column  (1 
pixel  from  a  block  of  16),  and  1  pixel  in  8  for  each  row  and  column  (1  pixel  from  a  block 
of  64  pixels).  The  sub-sets  were  used  to  assess  the  accuracy  of  data  reconstruction  by  the 
two  methods.  Table  5  gives  the  summary  statistics  for  the  full  data  set. 


Table  5:  Summary  statistics  for  NIR  for  the  128  pixels  by  128  pixels  region  of  Fort  A.  P. 
Hill 
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Figure  6;  Pixel  map  of  the  near  infra  read  (NIR)  of  part  of  the  SPOT  image  (128  by 
128  pixels)  for  Fort  A.  P.  Hill. 
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Wavelet  analysis 

The  method  of  wavelet  analysis  that  Edward  Bosch  used  was  that  of  Daubechies 
wavelets  (Daubechies,  1988).  The  size  of  the  image  was  f  x  ,  i.e.  128  rows  and  128 
columns  of  pixel  information  for  the  NIR  waveband.  This  region  chosen  from  the  part 
of  the  SPOT  image  that  we  analysed  and  described  in  the  previous  report  (1998),  and 
the  size  was  such  to  avoid  any  need  to  pad  the  data  to  create  appropriate  resolution 
levels  for  the  wavelet  analysis. 

The  wavelet  transform  was  done  with  a  pair  of  filters  and  not  the  wavelets  themselves. 
The  two-dimensional  wavelet  transform  is  separable  which  means  that  we  can  apply 
the  filters  in  the  horizontal  direction  and  then  in  the  vertical  direction  to  obtain  the 
desired  results.  As  we  described  above,  the  convolutions  can  be  represented  in  terms 
of  two  matrices  L  and  H.  The  matrix  L  is  made  up  of  shifts  of  the  low  frequency  filter 
l[n],  and  H  is  made  up  of  shifts  of  the  high  frequency  filter  h[n].  These  shifts  force 
every  other  point  in  the  output  of  each  convolution  to  be  discarded  (decimated). 

In  this  study  we  used  Daubechies  6-component  wavelet  filters.  This  wavelet  satisfies 
three  orthogonality  conditions  and  three  vanishing  moments.  The  orthogonality 
condition  implies  that  the  filter  is  orthogonal  to  its  shifts  by  two.  Increasing  the 
number  of  vanishing  moments  of  the  wavelet  and  its  filter,  forces  the  wavelet  to  be 
smoother  (continuous  and  differentiable),  but  at  the  same  time  it  increases  the  support 
of  the  wavelet  and  the  number  of  sampling  points  in  the  corresponding  filters. 

As  we  mentioned  above,  the  wavelet-coefficient  yvj{s,u)  provides  a  measure  of 
variation  of  the  signal  f(t)  at  s  and  u.  The  smoother  the  signal,  with  three  vanishing 
moments,  the  high  pass  filter  zeros  out  the  low  frequency  content  in  the  signal  leaving 
behind  little  information.  Since  the  high  pass  filter  h[n]  is  designed  to  discard  a  certain 
amount  of  information,  which  depends  on  the  number  of  vanishing  moments,  the  low 
pass  filter  l[n]  thus  must  retain  most  of  the  ‘energy’  or  information. 

Given  that  the  image  used  is  128  rows  by  128  columns  in  size  x  2’),  using  Mallat’s 
scheme,  the  original  resolution  level  is  0.  Applying  the  discrete  wavelet  transform 
(DWT)  to  the  data  once  produces  high  and  low  frequency  components  at  resolution 
level  1.  This  results  in  four  quarter  sets  of  data  of  size  64  by  64  pixels.  The  first 
quarter  of  the  data  represents  in  essence  a  sample  of  1  in  2  of  the  rows  and  columns  of 
the  data  matrix.  This  quarter  contains  the  scaling  coefficients  which  correspond  to  the 
low  pass  filter.  The  other  three  quarters  contain  the  wavelet  coefficients  which  are  high 
frequency;  quarter  two  contains  the  vertical  coefficients,  quarter  3  the  horizontal  ones 
and  quarter  4  the  diagonal  coefficients. 

For  this  analysis  the  low  and  high  frequency  components  were  used  separately  to 
assess  their  overall  contribution  to  the  original  image  and  to  analyse  the  effect  this  had 
on  the  variogram  of  the  original  image.  That  is,  we  computed  the  inverse  wavelet 
transform  at  several  levels  with  only  the  corresponding  low  frequency  components. 
Also,  we  computed  the  inverse  wavelet  transform  at  several  levels  without  the  low 
frequency  content. 
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Geostatistical  analyses 

The  variogram  was  computed  and  modelled  as  usual  using  the  smaller  data  set  with  a 
total  of  16  384  pixels,  Figure  7.  This  variogram  was  then  used  with  the  pixel 
information  to  filter  the  information  by  factorial  kriging  into  the  long-range  and^ort- 
range  components.  Ordinary  kriging  was  used  to  estimate  the  values  of  MR  at 
positions  where  pixels  had  been  removed  from  the  data.  In  other  words  the  estirnates 
coincided  with  the  locations  of  the  original  values  so  that  a  direct  comparison  could  be 
made  between  the  estimates  and  these. 


Results 

The  variogram  for  the  new  data  is  still  a  nested  structure,  but  the  correlation  ranges  are 
smaller  than  for  the  larger  part  of  the  scene  that  we  investigated  before.  The  model 
fitted  was  a  nested  spherical  function  with  two  structures.  Since  the  variogram  was 
somewhat  wavy  at  the  longer  lags,  to  improve  the  fit  I  modelled  it  to  a  lag  of  40  on  y. 
The  short-range  structure  was  6.6  pixels  or  130  m  and  the  long-range  structure  was 
pixels  or  420  m.  The  experimental  variogram  (points)  and  the  fitted  model  (line)  are 
given  in  Figure  7  a.  The  parameters  of  the  models  fitted  to  100  and  to  40  lags  are  given 
in  Table  6.  Since  the  data  were  skewed  I  transformed  them  using  Hermite  polynornials 
and  computed  the  variogram  from  the  transformed  data.  Figure  7  b  and  Table  6.  There 
was  little  difference  from  the  raw  variogram,  therefore,  I  did  the  analyses  on  the  raw 

data. 
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Multiresolution  analysis  (filtering) 


The  variogram  suggests  that  there  are  two  clear  scales  of  spatial  variation  present:  one 
of  about  120  m  and  the  other  of  about  420.  This  is  also  evident  in  the  pixel  map  of  the 
ordinary  kriged  estimates.  Figure  8.  There  is  local  detailed  variation  superimposed  on  a 
broader  pattern  of  variation.  The  major  large  structures  in  the  variation  that  are  evident 
appear  to  be  related  to  major  relief  forms:  the  drainage  basins  and  the  intervening 
spurs,  and  the  major  types  of  ground  cover.  Short-range  variation  is  also  evident 
related  to  the  water  bodies,  buildings  and  the  more  local  changes  in  ground  cover  and 
drainage.  These  were  described  in  the  previous  final  report. 


Table  6:  Model  parameters  for  the  variograms  computed  for  the  128  by  128  pixel  area 
of  the  SPOT  image 


Variable 

Model  type 

Nugget 

variance 

Sill(l) 

variance 

Range(l) 
pixels  (m) 

Sill(2) 

variance 

Range(2) 

pixels 

(m) 

MR 

(100  lags) 

Nested 

Exponentia 

1 

0 

227.0 

13.2 

(264) 

4.40 

85.1 

(1701.6) 

MR 
(40  lass) 

Nested 

Spherical 

0 

152.2 

nIEni 

91.71 

IBH 

Nested 

Spherical 

0 

0.5240 

mmlM 

0.3755 

IBH 

Long-range 

comnonent 

Circular 

0 

113.9 

ism 

0 

87.9 

4.30 

(86) 

0 

146.7 

8.40 

(168) 

74.82 

mSSM 

High  frequ 
-ency  1  in  2 

Pure  nugget 

High  frequ 
-ency  1  in  4 

Circular 

0 

41.7 

2.944 

(58.88) 

High  frequ 
-ency  1  in  8 

Circular 

0 

68.30 

6.690 

(133.8) 
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Figure  8:  Pixel  map  of  the  kriged  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels) 
for  Fort  A.  P.  Hill 
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Factorial  kriging  and  wavelet  analysis  enable  the  diflferent  spatial  scales  to  be  separated  in 
theory.  For  factorial  kriging  this  is  controlled  by  the  variogram  which  describes  the  variation 
present  in  the  data.  For  the  wavelet  analysis  this  is  controlled  by  the  resolution  level,  which 
in  turn  is  controlled  by  the  octave  bands  and  as  a  consequence  is  more  arbitrary.  One  aspect 
of  future  research  is  to  consider  how  the  variogram  could  be  used  to  guide  the  factoring 
process  that  controls  the  wavelet  multiresolution  analysis. 


Figure  9  is  a  pixel  map  of  the  kriged  estimates  of  the  long-range  component  of  the  variation 
filtered  using  the  variogram.  The  large  scale  variation  is  related  to  the  main  relief  features. 
The  band  of  dark  colours  in  the  North  and  central  part  of  the  map  are  damper  and  wetter 
areas,  and  the  lighter  ones  the  spurs,  upper  slopes  and  built  areas.  This  map  could  be  used 
effectively  to  guide  future  sampling.  If  the  end-user  is  interested  in  retrieving  this  level  of 
information  then  a  suitable  sampling  interval  can  be  chosen  using  the  range  of  the 
variogram.  A  sampling  interval  of  200  m  would  be  adequate  to  ensure  that  this  resolution 
of  variation  is  identified. 

The  pixel  map  of  the  short-range  variation  (Figure  10)  shows  the  detail  that  is  also  evident 
in  Figure  8,  but  less  clearly  so.  The  lakes  are  recovered  well  by  this  resolution.  The  dark 
patches  in  the  bottom  left  hand  comer  (1  to  20  on  the  x-axis  and  55  to  60  on  the  y-axis),  in 
the  central  area  (64  to  90  on  the  x-axis  and  1 15  to  125  on  the  y-axis),  and  at  the  top  of  the 
map  (45  to  70  on  the  x-axis  and  158  to  180  on  the  y-axis).  The  road  running  N-S  is  also 
evident  extending  N  along  longitude  100  (on  this  map).  The  other  short-range  stmctures 
probably  relate  to  changes  in  local  drainage  conditions  and  vegetation.  For  many  surveys 
recovering  this  intensity  of  variation  at  a  scale  of  about  120  m  would  require  too  much 
sampling.  A  sampling  interval  of  50  m  to  60  m  would  be  needed  to  resolve  this  short-range 
variation.  If  a  sampling  scheme  of  about  200  m  were  recommended  in  relation  to  the  long- 
range  variation  this  information  on  short-range  variation  would  be  lost.  These  maps  enable 
us  to  demonstrate  to  the  end-user  the  extent  of  information  that  is  likely  to  be  lost  by 
adopting  the  coarser  sampling.  Sampling  between  60  m  and  200  m  would  be  of  little  benefit 
because  most  of  the  short-range  variation  would  not  be  identified  and  sampling  at  less  than 
200  m  would  be  inefficient  to  identify  the  long-range  variation.  Variograms  were  computed 
fi'om  the  estimates  of  the  long-range  (Figure  11a)  and  the  short-range  (Figure  11  b) 
components.  They  recover  the  spatial  scale  of  the  variation  quite  well,  but  both  variograms 
were  difficult  to  model  satisfactorily. 

For  the  first  wavelet  analysis  the  level  of  resolution  was  1 .  The  coefficients  were  derived  as 
described  earlier.  The  low  frequency  and  the  high  frequency  coefficients  were  reconstructed 
by  the  inverse  wavelet  transform,  which  restored  each  of  the  64  by  64  sets  coefficients  to 
the  size  of  the  original  data  set.  These  are  shown  as  pixel  maps  and  should  be  compared 
with  the  appropriate  kriged  and  the  filtered  maps.  Figures  8  to  10.  In  addition  variograms 
were  computed  for  each  of  these  four  reconstructions:  one  low  frequency  (Figure  12  a)  and 
three  high  frequency  ones  (the  one  is  shown  in  Figure  12  b  is  the  average  as  they  were 
similar). 

The  low  frequency  reconstruction.  Figure  13,  is  veiy  similar  to  the  ordinary  kriged  output 
for  the  image.  Figure  8.  It  is  important  to  remember  that  the  ordinary  kriged 
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Figure  9;  Pixel  map  of  the  long-range  component  of  the  variation  in  NIR  of  part  of  the 
SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Short-range  estimates  of  NIR 
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Figure  10:  Pixel  map  of  the  short-range  component  of  the  variation  in  NIR  of  part  of 
the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  11:  a)  Variogram  of  the  long-range  component,  and  b)  variogram  of  the  short- 
range  component  of  the  variation  in  NIR  of  part  of  the  SPOT  image  (128  by  128 
pixels)  for  Fort  A.  P.  Hill 
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map  was  made  from  estimates  using  all  of  the  data,  whereas  the  low  frequency  wavelet 
reconstruction  used  the  1  in  2  sample,  i.e.  25%  of  the  original  data.  Both  the  long-  and 
short-range  components  of  the  variation  are  evident,  although  there  has  been  some  loss 
of  detail  in  the  short-range  variation.  For  example  the  road  is  less  clear  in  Figure  13 
than  in  Figure  8.  The  variogram  computed  from  the  low  frequency  reconstruction. 
Figure  12  a,  was  very  similar  to  the  variogram  of  the  raw  data.  Figure  7  a.  Hence  the 
spatial  structure  at  both  scales  has  been  retained  at  this  level  of  resolution.  The  most 
surprising  finding  was  that  related  to  the  high  frequency  reconstruction.  The  map  of 
the  high  frequency  component  (Figure  14)  does  not  appear  to  reflect  the  kind  of 
variation  present  in  the  map  of  the  short-range  component  from  factorial  kriging. 
However,  when  we  examined  them  in  detail  there  is  some  weak  evidence  of  the  lakes, 
which  are  so  clear  in  Figure  10.  The  variograms  computed  from  these  data  are  pure 
nugget,  as  shown  in  Figure  12  b.  This  means  that  the  high  frequency  components  are 
noise  at  this  level  of  resolution;  they  contain  no  spatial  structure.  The  latter  is  all 
retained  in  the  low  frequency  reconstruction. 

To  determine  whether  we  could  retrieve  the  long-  and  short-range  components  using 
wavelets  we  explored  the  next  resolution,  2,  in  effect  a  sampling  of  1  in  4  (or  1  pixel  in 
16).  Figure  15  shows  the  low  frequency  reconstruction.  There  is  still  long-  and  short- 
range  variation  evident,  although  the  short  range  variation  is  becoming  less  distinct;  for 
example  the  road  and  the  lakes  are  still  visible  but  their  margins  are  less  clearly 
defined.  Figure  16  shows  the  pixel  map  for  the  average  of  the  high  frequency 
reconstruction  and  it  is  clear  that  there  is  more  of  the  short-range  component  of  the 
variation  evident.  The  variogram  of  the  high  frequency  reconstruction  now  shows 
some  structure.  Figure  17  a.  Table  6  gives  the  model  parameters  of  this  variogram. 

The  low  frequency  reconstruction  of  the  1  in  8  resolution  3  now  shows  the  long-range 
component  of  the  variation  identified  by  factorial  kriging.  Figure  18.  This  resolution  is 
fairly  close  to  the  short-range  component  of  the  variogram,  i.e.  6.5  pixels,  and  this 
level  of  variation  appears  to  have  been  filtered  out  now.  So  it  seems  that  once  the 
resolution  of  the  short-range  structure  has  been  reached  the  effect  was  to  remove  the 
short-range  variation.  The  map.  Figure  19,  of  the  high  frequency  reconstruction  now 
shows  some  of  the  features  evident  in  the  kriged  map  of  the  short-range  component  of 
the  variation.  In  particular  the  lakes  are  evident.  The  variogram  computed  from  the 
average  of  the  three  high  frequency  reconstructions.  Figure  17  b,  shows  clear  evidence 
of  structure  and  the  range  of  spatial  correlation  described,  6.69  pixels,  is  close  to  the 
short-range  component  of  the  variogram  of  the  original  data. 


Summary 

It  is  clear  that  factorial  kriging  works  well  with  multiresolution  data.  The  main  reason 
for  this  is  that  the  filtering  is  controlled  by  the  variogram  which  is  a  function  of  the 
data  being  analysed.  It  is  a  valuable  method  for  directing  future  sampling  for  ground 
surveys  because  it  can  show  what  degree  of  variation  is  likely  to  be  recovered.  The 
multiresolution  analysis  using  wavelets  produces  a  different  outcome.  At  the  first 
resolution  the  high  frequency  components  remove  the  noise,  i.e  spatially  uncorrelated 
variation,  but  none  of  the  short-range  variation  is  resolved.  At  the  subsequent 
resolutions  close  to  the  short-range  component  in  the  variogram  the  high  frequency 
wavelets  recovered  the  short-range  variation,  especially  at  level  2,  evident  in  Figure  10. 
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10.  To  recover  the  long-range  component  of  variation  it  seems  that  choosing  a 
resolution  near  to  that  of  the  short-range  component  identified  by  the  variogram  is  an 
effective  way  of  avoiding  several  iterations  of  the  wavelet  analysis.  Once  the  scales  of 
variation  have  been  identified  by  the  variogram  the  choice  of  wavelet  coefficients 
retained  can  be  optimized.  This  presents  an  interesting  consideration  for  further  work. 
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Figure  13:  Pixel  map  of  the  low  frequency  reconstmction  from  the  wavelet  analysis  of 
NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill  at  a  resolution 
of  1  in  2 
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Figure  14:  Pixel  map  of  the  average  reconstruction  of  the  high  frequency  wavelets  ■ 

from  the  wavelet  analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for 

Fort  A.  P.  Hill  at  a  resolution  of  1  in  2  ^ 
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Wavelet  reconstruction  for  1  in  4  selection 
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Figure  15:  Pixel  map  of  the  low  jfrequency  reconstniction  from  the  wavelet  analysis  of 
MIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill  at  a  resolution 
of  1  in  4 
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Figure  16:  Pixel  map  of  the  average  high  frequency  reconstruction  from  the  wavelet 
analysis.of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill  at  a 
resolution  of  1  in  4 
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Figure  17:  a)  Variogram  of  the  high-frequency  component  for  the  1  in  4  resolution, 
and  b)  variogram  of  the  high-frequency  component  for  the  1  in  8  from  the  wavelet 
analysis  of  NIR  at  Fort  A.  P.  Hill 


40 


ABOVE  135 
130-135 
125-130 
120-125 
115-120 
110-115 
105-110 
100-105 
95-100 
90-  95 
85-  90 
BELOW  85 


Wavelet  reconstruction  for  1  in  8  selection 
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Figure  18:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet  analysis  of 
NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill  at  a  resolution 
of  1  in  8 
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High  frequency  reconstruction  1  in  8 
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Figure  19:  Pixel  map  of  the  average  high  frequency  reconstruction  from  the  wavelet 
analysis  of  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill  at  a 
resolution  of  1  in  8 
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Data  reconstruction 

The  128  by  128  pixels  were  sampled  by  taking  one  pixel  in  every  two  for  each  row  and 
column  (or  a  sample  of  1  in  4)  which  matches  resolution  level  1  in  the  wavelet 
analysis,  one  pixel  in  every  four  for  every  row  and  column  (or  a  sample  of  1  in  16),  a 
wavelet  resolution  of  2,  and  one  pixel  in  every  eight  for  each  row  and  column  (or  a 
sample  of  1  in  64),  a  wavelet  resolution  of  3.  The  low  frequency  wavelet  coefficients 
were  inverted  to  reconstruct  the  image  as  before.  Kriged  estimates  were  made  to 
coincide  with  the  original  data  points  for  each  data  set  using  the  variogram  model  from 
the  full  set  of  data.  These  maps  are  all  shown  as  pixel  maps. 

To  evaluate  the  accuracy  of  the  estimates  by  the  wavelet  reconstruction  and  kriging 
every  value  was  compared  with  the  original  values  of  NIR.  First  the  differences  were 
calculated  between  the  estimates  and  original  values  for  both  analyses  and  for  the  three 
sub-samples,  and  these  are  shown  as  pixel  maps  (Figures  21  and  22,  26  and  27, 31  and 
32).  The  statistical  distribution  of  these  differences  or  errors  has  also  been  determined 
and  these  are  shown  as  histograms  (Figures  23,  28  and  33).  In  addition  the  mean 
squared  differences  or  mean  squared  error  (mse)  was  calculated  (Table  4). 


Results 

The  results  were  not  entirely  what  we  expected  and  we  have  been  making  sure  that  the 
kriging  program  and  analyses  have  been  correct.  From  the  theory  of  geostatistics  we 
should  expect  that  the  kriged  estimates  would  have  the  smallest  mse,  but  they  do  not 
for  any  of  the  analyses.  It  was  this  that  led  us  to  explore  the  differences  in  more  detail 
to  try  to  gain  insight  into  the  results  from  the  two  methods.  In  spite  of  the  fact  that  the 
forward  and  inverse  wavelet  transform  are  linear  operators,  reconstructing  the  data 
with  only  some  of  the  wavelet  transform  components  is  not  done  in  a  linear  setting. 
Mallat  (1998)  says,  “It  is  often  easy  to  find  a  basis  that  produces  a  smaller  non-linear 
error  than  a  Karhunen-Loeve  basis,  ...”.  Although  in  this  study  we  are  not  using  the 
Karhunen-Loeve  transform  (principal  component  analysis),  some  methods  will  provide 
smaller  errors  than  others  depending  on  the  model  used. 


Sample  of  1  in  2 

The  pixel  maps  for  the  low  frequency  wavelet  reconstruction  and  kriging  from  the  1  in 
2  data.  Figures  13  and  20,  respectively  appear  to  be  very  similar  to  each  other.  The 
slight  ‘spottiness’  evident  on  the  kriged  map  is  because  punctual  kriging  was  used  and 
this  is  a  true  estimator  returning  the  value  at  the  data  points.  Table  7  gives  the  mean 
squared  errors  for  both  methods.  That  for  the  wavelets  is  less.  The  maps  of  errors  or 
comparisons.  Figures  21  (wavelet)  and  22  (kriging),  show  a  similar  pattern  in  general. 
However,  the  differences  between  them  help  to  explain  why  the  mse  is  greater  for  the 
kriged  estimates  than  for  the  wavelet  reconstruction.  There  eire  large  differences 
associated  with  the  lakes  where  there  are  clearly  marked  local  changes  associated  with 
boundaries  in  the  variation.  This  is  evidence  of  local  non-stationarity  which  violates 
the  assumptions  of  kriging.  Wavelets  are  known  to  be  suitable  for  dealing  with  local 
non-stationarity,  and  these  results  support  this.  Kriging  has  the  largest  absolute 
differences  and  there  are  more  of  them  than  for  the  low  frequency  wavelet 
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Figure  20:  Pixel  map  of  the  kriged  estimates  for  the  1  in  2  sample  of  NIR  of  part  of  the 
SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Comparisons  for  kriged  estimates  for  1  in  2 
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Figure  21:  Pixel  map  of  the  comparisons  between  the  kriged  estimates  for  the  1  in  2 
data  with  the  original  NIR  values  of  part  of  the  SPOT  image  (128  by  128  pixels)  for 
Fort  A.  P.  Hill 
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Comparisons  for  wavelet  reconstruction  1  in  2 


Figure  22:  Pixel  map  of  the  comparisons  between  the  low  frequency  wavelet 
reconstructed  values  for  the  1  in  2  data  with  the  original  NIR  values  of  part  of  the 
SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  23:  Histograms  of  a)  the  kriged  errors  and  b)  the  wavelet  errors  for  the  1  in  2 
sampling  for  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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reconstruction.  However,  compared  to  the  number  of  pixels  in  the  data  these  larger 
differences  are  few  compared  with  the  many  much  smaller  differences  for  the  majority 
of  the  estimates. 

To  explore  the  reasons  for  the  results  in  more  detail  the  histograms  of  the  differences 
were  examined.  Figure  23  a  and  b  are  the  histograms  of  the  differences  for  the  wavelet 
and  kriging  analyses,  respectively.  It  is  clear  that  punctual  kriging,  which  is  a  true 
estimator  at  the  data  points  has  a  larger  number  of  small  errors  than  the  wavelet 
analysis.  However,  this  is  not  consistent  as  the  number  of  data  points  retained  is 
reduced. 


Sample  of  1  in  4 

Figure  24  shows  the  result  of  kriging  this  sub-sample.  It  is  evident  that  much  of  the 
short-range  variation  has  been  lost  even  though  the  variogram  of  the  full  data  set  was 
used.  This  map  is  similar  to  that  for  the  long-range  component.  The  map  of  the  low 
frequency  reconstruction.  Figure  25,  shows  more  of  the  short-range  variation  and 
appears  to  be  much  more  accurate  visually  than  the  kriged  map.  The  maps  of  the 
differences.  Figures  26  and  27,  appear  to  be  similar  overall,  but  closer  examination 
shows  that  the  patches  where  the  differences  are  large  for  kriging  are  more  extensive 
than  those  for  the  wavelet  analysis.  The  values  of  the  MSEs  for  kriging  and  the 
wavelet  analysis  in  Table  7  also  suggest  that  kriging  performs  worse  than  wavelets. 
The  histograms.  Figure  28  a  (kriging)  and  b  (wavelets),  suggest  that  more  of  the 
kriged  values  have  smaller  differences  from  the  original  values  than  those  for  the 
wavelet  reconstruction.  However,  the  number  of  large  errors  is  also  greater  for  the 
kriged  values. 


Sample  of  1  in  8 

Figure  29  shows  the  result  of  kriging  this  sub-sample.  It  is  evident  that  much  more  of 
the  detail  in  the  variation  has  been  lost.  The  pattern  that  is  returned  is  coarse  and  no 
longer  reflects  even  the  long-range  component  of  the  variation  as  accurately.  Figure  30 
for  the  low  frequency  reconstruction  from  the  wavelet  analysis  also  shows  how  the 
detail  has  been  lost.  The  maps  of  the  differences.  Figures  31  and  32  are  again  similar, 
but  as  before  where  the  differences  are  greatest  for  the  kriged  differences  (Figure  30) 
so  their  extent  is  also  more  extensive.  It  seems  from  the  MSEs  Table  7  that  as  the  data 
become  more  sparse  and  separated  by  greater  distances  that  kriging  loses  power  in 
comparison  with  the  wavelet  analysis.  The  histograms.  Figure  33  a  (kriging)  and  b 
(wavelets),  suggest  that  the  wavelet  analysis  has  performed  better  with  this  sub-set  of 
the  data.  In  the  central  part  of  the  distribution  there  is  little  difference  between  the 
differences  for  wavelet  analysis  and  kriging,  but  there  seem  to  be  many  more  large 
errors  for  kriging  than  the  wavelets. 
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Figure  24:  Pixel  map  of  the  kriged  estimates  for  the  1  in  4  sample  of  NIR  of  part  of  the 
SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Wavelet  reconstruction  for  1  in  4  selection 
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Figure  25:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet  analysis  of 
NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill  at  a  resolution 
of  1  in  4 
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Figure  26:  Pixel  map  of  the  comparisons  between  the  kxiged  estimates  for  the  1  in  4 
data  with  the  original  NIR  values  of  part  of  the  SPOT  image  (128  by  128  pixels)  for 
Fort  A.  P.  Hill 
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Figure  27:  Pixel  map  of  the  comparisons  between  the  low  frequency  wavelet 
reconstructed  values  for  the  1  in  4  data  with  the  original  NIR  values  of  part  of  the 
SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  28:  Histograms  of  a)  the  kriged  errors  and  b)  the  wavelet  errors  for  the  1  in  4 
sampling  for  NCR.  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Kriged  estimates  for  1  in  8  selection 


25  50  75  100  125 


Figure  29:  Pixel  map  of  the  kriged  estimates  for  the  1  in  8  sample  of  NIR  of  part  of  the 
SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  30:  Pixel  map  of  the  low  frequency  reconstruction  from  the  wavelet  analysis  of 
NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill  at  a  resolution 
of  1  in  8 
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Figure  31:  Pixel  map  of  the  comparisons  between  the  kriged  estimates  for  the  1  in  8 
data  with  the  original  NIR  values  of  part  of  the  SPOT  image  (128  by  128  pixels)  for 
Fort  A.  P.  Hill 


Lv^iJ 


56 


I 


ABOVE 

20 

15- 

20 

10- 

15 

5- 

10 

0- 

5 

-5- 

0 

-10- 

-5 

-15- 

-10 

-20- 

-15 

BELOW 

-20 

Figure  32:  Pixel  map  of  the  comparisons  between  the  low  frequency  wavelet 
reconstructed  values  for  the  1  in  8  data  with  the  original  NIR  values  of  part  of  the 
SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Figure  33:  Histograms  of  a)  the  kriged  errors  and  b)  the  wavelet  errors  for  the  1  in  8 
sampling  for  NIR  of  part  of  the  SPOT  image  (128  by  128  pixels)  for  Fort  A.  P.  Hill 
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Summary 

The  histograms  of  the  differences  are  perhaps  the  most  illuminating  part  of  this 
analysis.  It  seems  that  we  need  to  explore  more,  but  that  kriging  performs  well  when 
fewer  data  have  been  removed  than  the  wavelet  analysis.  It  also  suggests  that  the  end- 
user  can  be  provided  with  some  insight  to  enable  them  to  choose  which  is  appropriate 
for  their  needs.  It  seems  that  for  the  1  in  2  and  lin  4  data  sub-sets  more  of  the  errors  are 
small  for  kriging  than  for  wavelets,  but  that  the  overall  error  is  least  for  the  wavelet 
analysis.  The  latter  is  clearly  more  successful  at  retmning  the  transition  features  present 
which  kriging  will  not  do  well.  Again  what  does  the  end  user  want? 

Another  thing  that  seems  to  emerge  from  this  investigation  is  that  the  variogram  could 
be  used  to  choose  an  optimal  subset  of  the  data,  based  on  the  distance  between  the 
values.  With  the  1  in  2  sample  both  the  long-range  and  short-range  components  of  the 
variation  were  restored  as  we  should  expect  from  the  correlation  structures  in  the 
variogram:  the  distance  between  the  pixels  was  less  than  the  range  of  the  short-range 
component.  With  the  1  in  4  sample  only  the  long-range  structure  is  successfully 
restored.  If  that  is  what  is  required  then  this  can  be  chosen  in  a  way  that  is  driven  by 
the  data  using  the  variogram. 

It  is  interesting  to  note  that  the  means  of  the  kriged  reconstructed  values.  Table  7,  are 
close  in  each  case  to  the  mean  of  the  original  data.  Table  5.  The  variances  for  the 
kriged  values  decrease  as  the  sampling  intensity  decreases  and  is  evidence  of  the 
smoothing  of  the  variation  that  occurs  with  kriging.  However,  the  variance  of  the 
kriged  values  for  the  1  in  2  sample  is  closer  to  the  original  variance  than  any  of  the 
other  analyses.  The  wavelet  analysis  retains  the  variance  better  as  the  sampling 
intensity  decreases.  Geostatistical  simulation  would  probably  perform  even  better  in 
terms  of  retaining  the  variance  in  the  data  and  this  method  should  also  be  compared 
with  wavelet  analysis  in  the  future. 
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Table  7;  Summary  values  from  comparisons  between  the  restored  values  and  the 
original  values  using  kriging  and  wavelets 


Data 

sub-sets 


Mean  difference 
(error) 

Mean 

squared  error 
(MSE) 

Mean  of 

reconstructed 

values 

Variance  of 

reconstructed 

values 

-0.0698 

28.361 

117.76 

235.64 

-0.3017 

87.122 

117.53 

196.97 

0.0627 

159.78 

117.89 

114.20 

-0.0959 

23.068 

122.77 

192.99 

0.0000 

67.860 

117.83 

210.13 

0.0000 

128.74 

117.83 

140.25 
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Part  IV:  The  analysis  of  the  vegetation  surveys  and  comparisons  with  the  SPOT 
data 

Introduction 

In  this  section  the  vegetation  surveys  and  analyses  will  be  described.  It  covers  the 
analysis  of  the  quantitative  data  from  surveys  1  and  2,  some  parts  of  these  data  sets 
have  been  combined,  and  the  qualitative  data  for  surveys  2  and  4.  The  analysis  of 
survey  3  was  included  in  the  final  report  for  the  previous  contract  (Contract  N68171- 
97-C-9029). 

Quantitative  Surveys  1  and  2 

Survey  1  was  carried  out  in  1997  at  A.  P.  Hill.  The  sample  comprises  several  small 
transects  that  have  random  starting  positions  within  the  seven  strata  of  the  training 
areas.  The  plot  size  corresponded  with  the  SPOT  pixel  size  of  20  m  by  20  m.  The 
points  along  the  transects  were  at  100  m  intervals  (see  Figure  34).  This  survey  mainly 
embraced  either  hard  or  soft  woodland  areas  of  vegetation.  The  second  survey  was  a 
square  grid  with  an  interval  of  300  m  covering  the  whole  of  our  study  site  at  A.  P.  Hill 
(Figure  35).  Since  there  were  many  sites  without  quantitative  woodland  information, 
because  it  included  grassland,  buildings  and  hard  standing,  the  sites  with  quantitative 
information  were  analysed  with  the  data  from  the  first  survey. 


Somple  1  (170) 


Figure  34:  Map  of  sites  for  Survey  1. 
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Sample  2  (114) 
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Figure  35:  Map  of  sites  for  Survey  2. 


Exploratory  data  analysis 

The  summaiy  statistics  of  the  17  quantitative  variables  were  analysed  for  surveys  1  and 
2  separately.  They  are  given  in  Tables  8  and  9.  The  skewness  values  are  generally  small 
showing  that  the  statistical  distribution  does  not  depart  seriously  from  normal,  except 
for  stem  spacing  (survey  1).  This  variable  had  one  extreme  value  which  was  removed  to 
obtain  a  near-normal  distribution  for  the  variogram  analysis.  Figures  36  and  37  show 
the  histograms  of  the  variables  listed  below  for  survey  1.  The  digital  numbers  for  the 
three  wavebands  of  the  SPOT  image  that  coincided  with  sites  where  the  vegetation  had 
been  examined  were  also  extracted  and  their  summary  statistics  are  given  in  Table  10 
for  both  surveys.  Their  histograms  are  shown  in  Figure  38. 

Variables  analysed  and  their  abbreviation: 

This  part  of  the  list  contains  those  variables  related  to  forest  density  (Set  A): 

maxcc  -  maximum  range  of  visual  estimate  of  crown  closure  (%) 
ovstmin  -  minimum  range  of  overstory  height  (ft) 
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ovstmax  -  mamimum  range  of  overstory  height  (ft) 
undstmn  -  minimum  range  of  understory  height  (ft) 
undstmx  -  maximum  range  of  understory  height  (ft) 
ba_f  -  estimate  of  basal  area  per  hectare  (metric  units) 
stem  *•  total  stems  in  plot  (count) 

ba_tot  -  sum  of  all  basal  area  for  each  tree  per  plot  (square  metres) 
stemsp  -  average  minimum  distance  between  stems  within  each  plot  (metres) 

This  part  of  the  list  contains  those  variables  related  to  tree  species  (Set  B): 

ba_so  -  percentage  of  total  basal  area  that  is  softwood  in  each  plot 
ba_ha  -  percentage  of  total  basal  area  that  is  hardwood  in  each  plot 
stem_so  -  percentage  of  total  number  of  stems  that  are  softwood  in  each  plot 
stem_ha  -  percentage  of  total  number  of  stems  that  are  hardwood  in  each  plot 
bad_so  -  percentage  of  dominant  basal  area  that  is  softwood  in  each  plot 
bad  ha  -  percentage  of  dominant  basal  area  that  is  hardwood  in  each  plot 
stemd_so  -  percentage  of  dominant  number  of  stems  that  are  softwood  in  each  plot 
stemd_ha  -  percentage  of  dominant  number  of  stems  that  are  hardwood  in  each  plot 


Table  8:  Summary  statistics  for  vegetation  measures  for  Survey  1 


Variable 

N 

Missing 

Mean 

Median 

Min 

Max 

Variance 

Standard 

deviation 

Skewness 

Kurtosis 

maxcc 

169 

1 

67.04 

70.0 

0.0 

100.0 

543.3 

23.31 

-1.22 

0.43 

minovst 

169 

1 

73.72 

80.0 

15.0 

110.0 

395.1 

19.88 

-1.19 

1.18 

maxovst 

169 

1 

78.54 

80.0 

20.0 

110.0 

405.1 

20.13 

-1.32 

1.29 

minunst 

169 

1 

11.35 

10.0 

0.0 

30.0 

33.3 

5.77 

0.96 

2.31 

maxunst 

169 

1 

20.66 

20.0 

0.0 

35.0 

62.2 

7.89 

-0.70 

0.28 

ba  f 

168 

1 

34.36 

34.3 

2.4 

16.1 

199.4 

14.12 

0.05 

0.31 

stem 

169 

1 

20.44 

19.0 

0.0 

81.0 

112.9 

10.62 

1.96 

7.01 

bat  tot 

169 

1 

1.07 

1.1 

0.0 

2.4 

0.2 

0.45 

0.01 

0.31 

stemsp 

168 

2 

2.33 

2.2 

0.9 

7.3 

0.6 

0.78 

2.00 

9.10 

ba_so 

168 

2 

46.16 

43.5 

0.0 

100.0 

1487.1 

38.56 

0.13 

-1.58 

ba_ha 

168 

2 

53.54 

54.5 

0.0 

100.0 

1487.1 

38.56 

0.13 

-1.58 

stem_so 

138 

2 

41.01 

33.3 

0.0 

100.0 

1338.9 

36.59 

0.33 

-1.41 

stem_ha 

168 

2 

58.99 

66.1 

0.0 

100.0 

1338.9 

36.59 

0.33 

-1.41 

bad_so 

168 

2 

48.94 

46.7 

0.0 

100.0 

1610.3 

40.13 

0.06 

-1.64 

bad_ha 

168 

2 

51.06 

54.3 

0.0 

100.0 

1610.3 

40.13 

0.06 

-1.64 

Stemd_ 

168 

2 

49.14 

48.8 

0.0 

100.0 

1570.7 

39.63 

0.02 

-1.62 

so 

Stemd_ 

ha 

168 

2 

50.86 

51.2 

0.0 

100.0 

1570.7 

39.63 

0.02 

-1.62 

Table  9:  Summary  statistics  for  vegetation  measures  for  Survey  2 


Variable 

N 

Missing 

Mean 

Median 

Min 

Max 

Variance 

Standard 

deviation 

Skewness 

Kurtosis 

maxcc 

60 

54 

68.17 

80.0 

5.0 

100.0 

674.5 

25.97 

-0.99 

-0.02 

minovst 

0 

114 

* 

* 

9)t 

* 

* 

« 

* 

* 

maxovst 

60 

54 

73.75 

80.0 

20.0 

100.0 

315.8 

17.77 

-1.32 

1.52 

minimst 

17 

97 

8.59 

10.0 

1.0 

20.0 

34.9 

5.91 

0.28 

-1.07 

maxunst 

54 

60 

15.11 

15.0 

3.0 

25.0 

23.9 

4.89 

-0.07 

-0.31 

ba  f 

58 

56 

32.06 

34.5 

3.4 

57.9 

175.9 

13.26 

-0.35 

-0.73 

stem 

58 

56 

19.91 

18.5 

5.0 

44.0 

96.1 

9.80 

0.63 

-0.29 

bat  tot 

58 

56 

1.07 

1.1 

0.1 

1.8 

0.17 

0.42 

-0.35 

-0.73 

stemsp 

58 

56 

2.50 

2.4 

1.2 

5.0 

0.57 

0.76 

1.05 

1.19 

ba__so 

58 

56 

57.54 

63.5 

0.0 

100.0 

1360.3 

36.88 

-0.31 

-1.46 

ba__ha 

58 

56 

42.46 

36.5 

0.0 

100.0 

1360.3 

36.88 

-0.31 

-1.46 

stem__so 

58 

56 

50.33 

52.1 

0.0 

100.0 

1288.4 

35.89 

-0.05 

-1.53 

stem_ha 

58 

56 

49.67 

47.9 

0.0 

100.0 

1288.4 

35.89 

-0.05 

-1.53 

bad_so 

58 

56 

60.87 

64.9 

0.0 

100.0 

1472.5 

38.37 

-0.38 

-1.43 

bad_ha 

58 

56 

39.14 

35.1 

0.0 

100.0 

1472.5 

38.37 

-0.38 

-1.43 

Stemd_ 

58 

56 

60.15 

71.8 

0.0 

100.0 

1520.1 

38.99 

-0.34 

-1.52 

so 

Stemd 

ha 

58 

56 

39.93 

28.2 

0.0 

100.0 

1520.1 

38.99 

-0.34 

-1.52 

Table  10:  Summary  statistics  for  the  three  wavebands  from  the  SPOT  data  for  Surveys 
1  and  2 


Variable 

N 

Missing 

Mean 

Median 

Min 

Max 

Variance 

Standard 

deviation 

Skewness 

Kurtosis 

Red  (1) 

116 

54 

61.94 

61.0 

58.0 

80.0 

14.5 

3.81 

2.49 

7.18 

Green 

116 

54 

36.33 

34.0 

32.0 

67.0 

33.56 

5.79 

3.17 

11.71 

(2) 

NIR(3) 

116 

54 

119.3 

121.0 

62..0 

148.0 

249.5 

15.79 

-0.55 

0.57 
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Figure  36:  Histograms  of  variables  in  Set  A  of  Survey  1. 


stem  so 


102030405060708090 

102030405060708090 

10  20  30  40  50  60  70  80  90 

cr 

Q 

o_ 

1 

CO 

o 

bod  ha 

stemd_so 

102030405060708090 


102030405060708090 


10  20  30  40  50  60  70  80  90 


stemd_hQ 


stemsp 


102030405060708090 


0,81.62.43.24.04.85.66.47.2 


Figure  37:  Histograms  of  variables  in  Set  B  of  Survey  1 . 
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band1  /  Sample  1  band2  /  Sample  1 


bands  /  Sample  1 


Figure  38:  Histograms  of  wavebands  1  (Red),  2  (Green),  3  (Nik)  and  NDVI  for  sites 
coinciding  with  Surveys  1  and  2. 


To  assess  which  of  these  variables  were  likely  to  represent  the  variation  the  data  most 
strongly  a  principal  components  analysis  was  done  on  the  correlation  matrix.  The  latter 
was  used  because  it  effectively  standardizes  the  data.  The  first  component  accounted 
for  53.7%  of  the  variation  and  the  second  18%.  The  variables  that  ‘loaded'  most 
heavily  on  the  first  component  were; 


ba  so,  ba_ha,  stem  so,  stem  so,  stem_ha,  bad  so,  bad_ha,  stemd_so  and  stem  ha. 
The  variables  that  ‘loaded'  most  heavily  on  the  second  component  were: 
maxcc,  ba_f,  stem,  baJ:ot  and  stemsp. 
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A  set  of  variables  that  is  considered  to  express  the  variation  and  summarise  it 
adequately  is: 

maxcc,  ba_f,  stem,  stemsp,  ovstmax,  undstmx  and  ba_so. 

These  are  based  on  the  distribution  of  the  variables  in  the  plane  of  PCI  and  PC2  (Figure 
39). 


PCA  -  Loadings.  Sample  1 


Figure  39:  Plot  of  variables  based  on  their  loadings  in  the  plane  of  PCI  and  PC2. 

Table  11  gives  the  correlations  for  the  vegetation  measures  and  the  DNs  of  the 
wavebands.  In  general  these  are  small  for  the  vegetation  measures  and  DNs.  Those  for 
NER.  are  the  largest  for  maxcc,  stem_so,  ste_ha,  stemd_so  and  stem_ha.  There  are  some 
strong  correlations  for  the  vegetation  measure  which  are  to  be  expected,  for  example 
ba_so  and  ba_ha  which  add  to  100%. 
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Table  11.  Correlations  for  the  vegetation  measures  and  the  three  wavebands  from  the 
SPOT  image. 


***  Correlation  matrix  *** 


bandl 

1.000 

band2 

0.960 

1.000 

band3 

-0.205 

-0.308 

1.000 

cc 

0.023 

0.017 

0.227 

1.000 

ovstmin 

0.036 

0.018 

0.178 

0.205 

1.000 

ovstmax 

0.030 

0.021 

0.198 

0.237 

0.956 

1.000 

undstmn 

-0.148 

-0.113 

0.076 

0.145 

0;196 

0.217 

1.000 

undstmx 

0.073 

0.080 

0.099 

0.260 

0.427 

0.487 

0.515 

ba_f 

0.179 

0.177 

0.087 

0.518 

0.557 

0.594 

0.141 

stem 

0-081 

0.098 

-0.021 

0.474 

-0.346 

-0.306 

0.021 

ba_tot 

0.178 

0.177 

0.087 

0.518 

0.557 

0.594 

0.141 

ba„so 

0.005 

-0.005 

-0.182 

-0.260 

-0.596 

-0.625 

-0.014 

ba_ha 

-0.005 

0.005 

0.182 

0.260 

0.596 

0.625 

0.014 

stem_so 

0.004 

0.010 

-0.224 

-0.230 

-0.676 

-0.687 

-0.053 

stem_ha 

-0.004 

-0.010 

0.224 

0.230 

0.676 

0.687 

0.053 

bad_so 

0.002 

-0.015 

-0.158 

-0.276 

-0.541 

-0.572 

0.007 

bad_ha 

-0.002 

0.015 

0.158 

0.276 

0.541 

0.572 

-0.007 

stemd_so 

0.000 

-0.006 

-0.211 

-0.274 

-0.547 

-0.568 

0.009 

stemd_ha 

0.000 

0.006 

0.211 

0.274 

0.547 

0.568 

-0.009 

stemsp 

-0.007 

-0.033 

-0.030 

-0.456 

0.182 

0.136 

-0.015 

bandl 

band2 

band3 

cc 

ovstmin 

ovstmax 

undstmn 

undstmx 

1.000 

ba_f 

0.396 

1.000 

stem 

-0.014 

0.260 

1.000 

ba_tot 

0.396 

1.000 

0.260 

1.000 

ba_so 

-0.146 

-0.385 

0.267 

-0.385 

1.000 

ba_ha 

0.146 

0.385 

-0.267 

0.385 

-1.000 

1.000 

stem_so 

-0.207 

-0.385 

0.309 

-0.385 

0.941 

-0.941 

1.000 

stem_ha 

0.207 

0.385 

-0.309 

0.385 

-0.941 

0.941 

-1.000 

bad_so 

-0.116 

-0.376 

0.243 

-0.376 

0.990 

-0.990 

0.894 

bad_ha 

0.116 

0.376 

-0.243 

0.376 

-0.990 

0.990 

-0.894 

stemd_so 

-0.121 

-0.349 

0.269 

-0.349 

0.972 

-0.972 

0.939 

stemd_ha 

0.121 

0.349 

-0.269 

0.349 

-0.972 

0.972 

-0.939 

stemsp 

-0.174 

-0.302 

-0.660 

-0.302 

-0.122 

0.122 

-0.189 

undstmx 

ba_f 

stem 

ba_tot 

ba_so 

ba_ha 

stem_so 

stem_ha 

bad„so 

bad_ha 

stemd_so 

stemd_ha 

stemsp 

1.000 

-0.894 

0.894 

-0.939 

0.939 

0.189 

1.000 

-1.000 

0.966 

-0.966 

-0.076 

1.000 

-0.966 

0.966 

0.076 

1.000 

-1.000 

-0.137 

1.000 

0.137 

1.000 

stem_ha 

bad_so 

bad_ha 

stemd_so 

stemd_ha 

stemsp 
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Variogram  analysis 

Experimental  variograms  were  computed  for  all  of  the  variables  listed  above  for  the 
combined  data  from  surveys  1  and  2.  Variograms  were  computed  in  four  directions  at 
the  outset,  but  the  number  of  sites  is  marginal  for  this.  For  set  A  variables  the  directions 
of  maximum  and  minimum  variation  are  not  consistent,  but  for  set  B  variables  the 
variation  in  direction  NNE  to  SSW  (o)  have  the  longest  range  of  spatial  dependence  and 
the  largest  sill  variances  and  those  at  right  angles  have  the  shortest  ranges  and  the 
smaller  sill  variances  (*)  (Figures  40  and  41). 


Figures  42  and  43  show  the  experimental  omnidirectional  variograms  for  the  two  sets  of 
variables  from  surveys  1  and  2.  Those  that  show  reasonable  spatial  structure  are: 
maxcc,  ovstmin,  ovstmax,  undstmn,  stem,  ba_so,  ba_ha,  stem  so,  stem  ha,  bad_so, 
bad_ha,  stemd_so  and  stemd_ha.  For  the  twin  variables,  such  as  ba_so  and  ba_ha  the 
variograms  are  identical  for  the  reasons  given  earlier.  The  following  variables  were 
modelled:  maxcc,  overstory  height  (derived  from  ovstmin  and  ovstmax),  understory 
height  (derived  from  undstmn  and  undstmx),  ba_f,  stem,  stem  spacing,  ba_so 
(equivalent  to  ba  ha  also),  stem  so,  bad  so  and  stemd_so.  In  addition  the  multivariate 
variogram  from  this  analysis  was  computed  and  modelled,  also  elevation,  and  the  three 
wavebands  and  NDVI.  They  are  shown  in  Figures  44  to  47. 

Table  12  gives  the  model  parameters  of  the  variables  modelled.  The  experimental 
variograms  of  many  of  the  properties  in  Table  12  are  somewhat  erratic.  This  could  be 
related  to  the  irregular  sampling  scheme.  However,  there  appears  to  be  some  evidence 
of  periodicity  in  several  variograms  with  wavelengths  of  between  500  m  and  700  m.  A 
previous  report  that  contained  transects  of  the  pixels  to  match  the  vegetation  ones  also 
showed  periodicity  in  the  DNs.  There  appears  to  be  some  relation  between  the  range  of 
spatial  dependence  of  elevation  and  several  of  the  vegetation  measures.  The  multivariate 
variogram  has  identified  a  short  range  component  of  variation  of  just  over  300  m  which 
matches  with  the  short  range  component  of  NIR.  The  variograms  of  the  vegetation 
classes  are  described  later  in  this  report.  The  models  fitted  to  directional  variograms  of 
ba_so  are  revealing:  the  variation  in  direction  135°  is  462  m  and  that  in  direction  45°  is 
1271  m.  This  suggests  that  the  different  ranges  might  reflect  some  anisotropy  in  the 
variation.  This  was  identified  in  the  image  data,  but  because  the  sill  heights  were 
different  this  signalled  zonal  anisotropy  which  cannot  be  corrected  simply.  It  suggests 
that  there  are  distinct  strata  present  and  this  is  evident  from  the  areas  with  different 
kinds  of  vegetation.  There  are  also  distinct  landscape  units  which  will  be  explored  in 
the  next  report. 
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Figure  40:  Directional  experimental  variograms  for  Set  A  variables  for  Surveys  1  and 
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Figure  42;  Omnidirectional  experimental  variograms  of  Set  A  variables  from  Surveys 
1  and  2, 
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Figure  43:  Omnidirectional  experimental  variograms  of  Set  B  variables  from  Surveys  1 
and  2. 
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Table  12.  Variogram  model  parameters  for  the  quantitative  information  from  Surveys  1 
and  2,  and  for  elevation  and  the  image  information. 


Variables 

Model  type 

Nugget 

Sill 

Sill 

Range 

Range 

variance 

Cl 

C2 

aj  (m) 

02  (m) 

Canopy  closure 

Circular 

379.1 

211.4 

1707.0 

Overstory  height 

Double 

spherical 

0 

275.9 

124.9 

243.0 

2034.0 

Understory 

height 

Circular 

19.2 

12.6 

1391.0 

Basal  are  (field) 

Spherical 

69.9 

126.5 

232.0 

Stem 

Pentaspherical 

21.4 

65.7 

380.0 

Stem  spacing 

Circular 

0.305 

0.193 

407.0 

ba_so.ha 

Double 

spherical 

0 

980.2 

563.5 

182.0 

1553.0 

ba_so/ha  (45°) 

Circular 

662.3 

1271.0 

1271.0 

ba_so/ha  (135°) 

Circular 

428.0 

841.4 

462.0 

stem_so/ha 

Spherical 

892.7 

838.3 

1428.0 

bad_so/ha 

Spherical 

909.1 

819.9 

1274.0 

stemd_so/ha 

Circular 

839.9 

869.7 

1432.0 

Multivariate 

variogram 

Circular 

6.05 

3.08 

309.9 

Elevation 

Circular 

93.05 

312.5 

1562.0 

Red  (1) 

Pentaspherical 

2.21 

17.89 

906.0 

Green  (2) 

Double 

spherical 

0 

20.6 

19.7 

386.0 

1047.0 

NIR(3) 

Circular 

99.24 

145.0 

673.6 

NDVI 

Double 

spherical 

0.0015 

0.00307 

0.00202 

666.8 

1261.0 

Ccncpy  closure 
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Figure  44:  Experimental  variograms  and  fitted  models  for  Surveys  1  and  2. 
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Figure  46:  Experimental  variograms  and  fitted  models  Surveys  1  and  2. 
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Cross  variograms 

The  theory  for  computing  cross  variograms  between  two  or  more  variables  is  given  at 
the  beginning  of  the  report.  Cross  variograms  were  computed  between  the  vegetation 
measures  and  the  DNs  from  the  three  SPOT  wavebands.  Those  selected  and  shown  in 
Figures  47  to  50  show  some  relation  between  the  variables.  For  band  1  (Red)  there  is  a 
negative  relation  between  maxcc,  unstmn  and  stem,  and  a  positive  relation  between 
stem  spacing  (Figure  47).  The  relations  with  the  other  variables  is  not  clear.  For  band  2 
(Green)  there  are  clear  negative  relations  with  maxcc  and  stem,  and  a  positive  relation 
with  stem  spacing  (Figure  48).  For  band  3  (NIR)  there  are  positive  relations  between 
maxcc,  ovstmax,  stem  and  ba_f,  and  a  negative  relation  with  ba_so  (Figure  49).  Cross 
variograms  with  elevation  are  give  in  Figure  50.  Overall  their  relations  with  the 
vegetation  measures  are  weak. 
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Figure  47:  Cross  experimental  variograms  between  band  1  (Red)  and  selected 
vegetation  measures  for  Surveys  1  and  2.. 
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Figure  48:  Cross  experimental  variograms  between  band  2  (Green)  and  selected 
vegetation  measures  for  Surveys  1  and  2.. 
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Figure  49:  Cross  experimental  variograms  between  band  3  (NIR)  and  selected 
vegetation  meaSutes  for  Surveys  1  and  2.. 
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Figure  50:  Cross  experimental  variograms  between  elevation  and  selected  vegetation 
measures  for  Surveys  1  and  2.. 
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The  analysis  of  the  qualitative  data  for  Surveys  2  and  4 

In  addition  to  the  measured  values  for  the  woodland  areas.  Surveys  2,  3  and  4  also 
described  the  vegetation  in  categories  or  classes.  The  distribution  of  the  sites  for  all  of 
the  ground  cover  surveys  is  given  in  Figure  51,  and  that  for  Survey  4  alone  in  Figure 
52.  The  sampling  scheme  for  Survey  4  was  along  a  series  of  transects  with  a  sampling 
interval  of  10  m.  This  survey  was  confined  to  a  smaller  area  than  the  transects  of 
Survey  3  which  had  a  50  m  sampling  interval. 

In  the  original  classification  of  these  data  there  were  several  classes  with  few  sites. 
Therefore,  as  for  Survey  3  described  in  the  previous  final  report  (Contract  N68171-97- 
C-9029),  the  number  of  classes  was  reduced.  The  original  and  new  classes  are  given  in 
Table  13.  For  Survey  2  there  were  114  sampling  points  and  for  Survey  4  there  were 
321  points.  Figure  53  shows  the  spatial  distribution  of  the  eight  classes  in  relation  to 
the  sampling  points  for  Survey  2  and  Figure  54  shows  their  distribution  for  Survey  4. 
Some  sites  had  a  mixture  of  ground  cover,  for  example  grass/low  bushes,  and  these 
were  eliminated  fi'om  the  analysis.  For  Survey  2,  19  sites  were  excluded  and  so  the 
variogram  analyses  were  carried  out  using  the  remaining  95  sites.  For  Survey  4,  71 
sites  were  excluded  leaving  250  sites  for  the  analysis.  Table  14  summarises  the 
contents  of  the  classes. 
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Figure  51:  Map  of  sites  for  the  four  ground  cover  surveys. 
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Table  13:  Reclassification  of  the  classes  of  ground  cover  for  sites  Surveys  2  and  4. 


Sample  2 

Manmade 

Building/Asphalt 

Buildings 

Campgmd-like 

Site 

Gravel  Parking 
Lx)t 

Powerline  Grass 
Private-Landfill 
Road/Grass 

Grass 

Grass 

Grass  Field 
Grass  Field 
Tall  Grass 

Oak 

Oak 

Hard  Mix 
Hard  Mix 
Oak  Mix 
Poplar 

Mixed_For. 

Mixed_For. 

Pine 

Pine 

Plantation 
Young  Pine 
Private-Pine 

Wetland 

Wetland 

Private-Wetland 

Wetland 

wAVaterways 

Excluded 

Grass/Cemetery 


Grass/Edge  of 
Mixed 


Grass/Houses 
Grass/Low  Bushes 
Grass/Rd/Bldg 
Grass/Road 
Lake 

Maple/Pine 

Oak/Holly 

Oak/Pine 

Pine/Maple 

Pine/Mix 

Pine/Plantation 

Marshy  Wetland 

Private-Mix 

Private-Mixed 


Sample  4 

Manmade  works 
Asphalt 
Asphalt/Grass 
Asphalt/Grass/Soi 
1 

Concrete 

roof-Asphalt 

Field 

Field 

Grass 

Grass 

Grass  (Tall) 

Tall  Grass 

Hardwood 
Hardwood 
Forest  (H) 

Mixed  Forest 
Pine/Hardwood 
Pine/HW 
Forest  (M) 

Pine 

Pine 

Forest  (P) 

Forest 

Forest 

Shrub/S  crubAV  e 
tland 

Shrub/ScrubAVetl 

and 

Excluded 
Edge  Field/Forest 
Edge/Field 
Field/Gravel  Road 
For  Dist  (M) 
For/Tall  Grass 


Forest  (Wet/Dry) 
Forest  Disturbed 


Forest/Grass 

Forest/Road 

Forest^Scrub 

Forest/Tall  Grass 

Grass/Asphalt 

Grass/Asphalt/Ce 

ment 

Grass/Conifer 

Grass/Dirt 

Grass/Forest 

Grass/HW 

Grass/Soil 

HW/Marsh 

Marsh/HW 

Pine/Brush 

Pine/Grass 

Pine/Hard/Shrub 

Pine/HW/Marsh 

Pine/Shrub 

Shrub/Scrub 

ShrubScrub/Forest 

Tall  Grass/Field 

Wet  Forest/HW 
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Table  14:  Summary  of  ground  cover  classes  for  sites  Surveys  2  and  4. 


I  Survey  2 

Survey  4 

Ground  Cover  Class 

Number  of  samples 

Ground  Cover  Class 

Number  of  samples 

Manmade 

11 

Manmade  works 

20 

Grass 

16 

Field 

5 

Grass  field 

3 

Grass 

71 

Oak 

3 

Hardwood 

24 

Hard  Mixed  Forest 

12 

Mixed  Forest 

32 

Mixed  Forest 

8 

Pine 

11 

Pine 

34 

Forest 

76 

Wetland 

S 

Shrub/ scrub/wetland 

11 

Excluded 

19 

Excluded 

71 

I  Total  sites  114 

Total  sites  321  { 

Sample  4  (321) 


Figure  52:  Map  of  sites  for  Survey  4. 
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Variogram  analysis  of  categorical  data  for  Survey  2 

Experimental  variograms  were  computed  for  each  of  the  eight  ground  cover  classes. 
Figure  55.  The  ones  that  show  clear  structure  are  for  the  categories:  manmade,  grass, 
oak,  pine  and  wetland.  Considering  the  small  sample  size  the  evident  structure  in  these 
suggests  that  there  are  distinct  areas  of  ground  cover  type  that  can  be  identified  with 
relatively  little  sampling  effort.  The  variogram  for  oak  appears  to  have  a  nested 
structure.  These  variograms  have  not  been  modelled  at  present,  but  this  could  be  done 
in  the  next  phase  of  work.  From  the  graphs  the  average  range  of  spatial  dependence  is 
about  600  m.  This  fits  in  well  with  the  long  range  component  of  the  variation  in  MR 
for  the  original  part  of  the  SPOT  image  analysed.  Table  15. 

The  experimental  multivariate  variogram  based  on  the  classes  was  computed.  This 
analysis  summarises  the  average  rate  of  change  from  one  class  of  ground  cover  to 
another.  Figure  56  shows  the  experimental  variogram  with  the  fitted  model.  The 
variogram  has  a  large  nugget  variance  because  there  were  relatively  few  sites  and  the 
sample  spacing  of  300  m  was  large.  The  variation  that  has  not  been  detected  in  the 
change  from  one  type  of  ground  cover  to  another  is  embraced  by  the  nugget  variance. 
The  variogram  model  was  a  single  structure  with  a  range  of  spatial  dependence  of  575 
m  (Table  15).  This  corresponds  closely  with  the  long-range  structure  identified  in  the 
nested  variogram  of  MR  computed  from  the  original  part  of  the  SPOT  image  that  was 
analysed.  This  was  slightly  larger  than  the  area  used  for  the  wavelet  analysis  (Part  II  of 
this  report).  The  variogram  of  this  waveband  was  fitted  by  a  nested  spherical  model 
with  a  short  range  component  of  120  m  and  a  long  range  component  of  542  m  (Table 
15).  It  is  evident  that  the  spectral  information  is  identifying  the  major  changes  in  the 
different  kinds  of  ground  cover. 


Table  15.  Variogram  model  parameters  for  the  wavebands  in  the  original  SPOT  image 
and  the  qualitative  information  from  Surveys  2  and  4. 


Variables 

Model  type 

Nugget 

Sill 

Sill 

Range 

Range 

variance 

Cl 

C2 

ai  (m) 

a2  (m) 

Red  -  average 

Double 

spherical 

0.4661 

8.080 

18.432 

165.56 

518.94 

Green  -  average 

Double 

spherical 

2.026 

16.568 

28.437 

224.32 

650.48 

MR-  average 

Double 

spherical 

0.0 

151.39 

118.87 

120.56 

542.22 

MDVI-  average 

Double 

spherical 

0.0 

0.0343 

0.0138 

139.74 

643.35 

Ground  cover  2 

Circular 

0.431 

0.283 

575.0 

Ground  cover  4 
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spherical 

0.112 
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Figure  53:  Distribution  of  sites  for  each  category  of  ground  cover  for  Survey  2. 
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Figure  54:  Distribution  of  sites  for  each  category  of  ground  cover  for  Survey  4. 
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Figure  56:  Experimental  multivariate  variogram  (symbols)  based  on  the  ground  cover 
classes  for  Survey  2  and  the  fitted  model  (solid  line). 


Cross  variograms  were  computed  as  before  between  the  ground  cover  classes  and  each 
of  the  three  wavebands.  Figures  57  to  59  show  the  cross  variograms  for  Survey  2. 
Those  for  wavebands  1  and  2  are  similar.  For  the  Red  waveband  (1)  manmade,  grass, 
oak,  hard  mixed  forest,  mixed  forest,  pine  and  wetland  show  evidence  of 
coregionalization  with  it.  The  strength  of  this  cannot  be  judged  at  present  until  they 
have  been  modelled.  The  coregionalization  for  the  Green  waveband  (2)  is  as  above, 
except  for  grass.  For  NIR  (3)  the  coregionalization  is  between  oak,  hard  mixed  forest, 
mixed  forest  and  wetland.  The  relations  for  the  two  grass  categories  and  pine  are  not 
strong. 
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Figure  59:  Cross  experimental  variograms  between  band  3  (NIR)  and  the  ground 
cover  classes  for  Survey  2. 
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Variogram  analysis  of  categorical  data  for  Survey  2 


Experimental  variograms  were  computed  for  each  of  the  eight  ground  cover  classes  to 
a  maximum  lag  of 400  m.  Figure  60  and  1000  m.  Figure  61 .  They  are  more  erratic  than 
the  ones  for  Survey  2,  which  is  surprising  considering  the  large  number  of  sampling 
points.  However,  the  behaviour  of  the  variograms  probably  reflects  the  fact  that  the 
categories  of  ground  cover  repeat  themselves  along  the  transects,  this  would  cause  the 
waviness  that  is  apparent  in  many  of  them.  The  variograms  all  show  clear  structure, 
especially  over  a  short  distance  of  about  100  m.  The  average  range  of  these 
variograms,  based  on  where  they  start  to  flatten,  is  about  145  m.  This  relates  closely  to 
the  short  range  component  of  the  variation  in  NIR  for  the  original  part  of  the  SPOT 
image  analysed.  Table  15.  These  variograms  could  be  modelled  in  the  next  phase  of 
the  work. 

The  experimental  multivariate  variogram  based  on  the  classes  was  computed  as  before. 
Figure  62.  The  fitted  model  is  the  solid  line  the  figure.  The  variogram  was  fitted  best 
by  a  nested  spherical  model  with  a  short-range  structure  of  67  m  and  long-range  one  of 
689  m.  However,  the  appearance  of  the  experimental  variogram  suggests  that  the 
short-  range  component  should  be  at  about  150m  to  175  m.  Several  attempts  have 
been  made  to  fit  a  better  model,  but  the  erratic  nature  of  the  experimental 
semivariances  has  prevented  this.  This  variogram  mirrors  the  form  of  that  for  NIR 
closely  in  spite  of  the  fact  that  the  model  parameters  are  different  for  the  short-range 
component. 

Cross  variograms  were  computed  as  before  between  the  ground  cover  classes  and  each 
of  the  three  wavebands.  Figures  63  to  65  show  the  cross  variograms  for  Survey  4. 
Those  for  the  Red  and  Green  wavebands  are  similar.  The  coregionalization  appears  to 
be  weaker  for  this  analysis,  but  this  could  be  related  to  the  complex  nature  of  the 
individual  variograms.  The  variables  that  show  the  strongest  relation  with  wavebands  1 
and  2  are  manmade,  field,  grass,  forest  and  wetland.  For  NIR  (waveband  3)  the 
strongest  relations  are  with  manmade,  hardwood,  mixed  forest,  pine,  forest  and 
wetland. 

Cross  variograms  were  also  computed  between  the  ground  cover  categories  and 
elevation.  The  strongest  relations  are  with  field,  mixed  forest,  forest  and  wetland.  In 
the  previous  final  report  we  commented  on  the  fact  that  the  patterns  in  the  long-range 
component  of  the  variation,  in  particular,  appeared  to  have  a  strong  relation  with  the 
physiography.  These  results  appear  to  confirm  this. 


Manmade  works 


Field 


Figure  60:  Experimental  variograms  of  the  ground  cover  classes  for  Survey  4 
computed  to  a  maximum  lag  of 400  m. 
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Figure  61:  Experimental  variograms  of  the  ground  cover  classes  for  Survey  4 
computed  to  a  maximum  lag  of  1 000  m. 
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Summary  of  vegetation  analysis 


The  results  of  the  ground  cover  analyses  suggest  that  there  are  distinct  spatial  patterns 
in  the  variation;  a  short-range  one  related  to  local  changes  in  cover  of  about  150  m 
extent  and  a  longer  range  one  of  about  550  m  in  extent.  They  relate  to  the  patterns 
observed  in  the  imagery,  in  particular  for  the  NIR  waveband. 


Multivariate  variogram  of  classes  in  Sample  4 


Figure  62;  Experimental  multivariate  variogram  (symbols)  based  on  the  ground  cover 
classes  for  Survey  4  and  the  fitted  model  (solid  line). 
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Figure  63:  Cross  experimental  variograms  between  band  1  (Red)  and  the  ground  cover 
classes  for  Survey  4. 


100 
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Figure  66:  Cross  experimental  variograms  between  elevation  and  the  ground  cover 
classes  for  Survey  4. 


101 


Digital  Elevation  Analysis 

Digital  elevation  information  for  A.  P.  Hill  was  provided  for  most  of  the  area  of  the 
SPOT  image  that  we  have  been  working  with.  Part  of  the  southwestern  corner  is 
missing  because  this  are  is  beyond  the  confines  of  the  base.  The  data  were  on  a  5  m 
grid.  Table  16  gives  the  summary  statistics  of  the  elevation  data.  Variograms  were 
computed  along  the  rows  and  columns  of  the  grid  and  the  average  variogram  was  also 
estimated.  Figure  67  shows  the  variograms.  They  all  have  a  concave  upward  slope  near 
to  the  origin  which  suggest  that  local  trend  or  drift  might  be  present.  This  is  common 
for  elevation  data. 


Table  16:  Summary  statistics  for  Elevation 


Elevation 

Number  of 

26878 

observations 

Minimum 

155.0 

Maximum 

232.0 

Mean 

198.15 

Variance 

234.88 

Standard  deviation 

15.33 

Skewness 

-0.5499 

To  assess  whether  the  local  trend  could  be  removed  by  a  coarser  sampling  interval  the 
data  were  sampled  to  produce  a  20  m  grid.  This  corresponded  with  the  information 
from  the  SPOT  image  with  a  pixel  size  of  20  m  by  20  m.  The  variograms  from  these 
data  now  show  no  evidence  of  local  trend,  but  there  is  evidence  of  regional  trend,  i.e. 
the  variograms  start  to  rise  after  an  initial  sill  has  been  reached  at  about  a  lag  distance  of 
800  m.  Figure  68.  The  amount  of  trend  present  was  determined  by  fitting  linear, 
quadratic  and  cubic  functions  to  the  coordinates  of  these  data.  The  linear  function 
accounted  for  17%  of  the  trend,  the  quadratic  22%  and  the  cubic  31%.  Variograms 
were  computed  from  the  residuals  to  these  trend  functions.  Figure  69  shows  the 
variograms  computed  from  the  residuals  of  the  linear  trend.  They  show  clearly  that  the 
regional  trend  has  been  removed.  The  variograms  computed  from  the  residuals  of  the 
quadratic  trend.  Figure  70,  are  similar,  as  are  those  for  the  cubic  residuals.  Figure  71. 
The  variogram  that  is  most  different  is  the  one  for  the  columns  in  Figure  71  (cubic). 
This  emphasises  the  possibility  of  some  periodicity  in  the  N-S  direction.  Table  17  gives 
the  parameters  of  the  models  fitted  to  the  average  variograms.  The  model  for  the  raw 
data  was  fitted  to  a  lag  of  600  m  only  because  thereon  the  variogram  continues  to 
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rise.  A  single  pentaspherical  model  provided  the  best  fit.  The  model  has  a  range  o 
578  m  which  relates  closely  to  the  long-range  component  in  the  image.  The  variograms 
of  the  residuals  were  all  fitted  best  by  nested  spherical  models.  Table  17.  Their  ranges 
do  not  correspond  as  closely  with  image  data  as  the  vegetation  does.  The  short-range 
component  of  about  220  m  is  similar  to  that  of  the  Green  waveband  (2)  and  the  long- 
range  one  of  about  1000  m  is  similar  to  that  fitted  to  the  rows  of  NDVI. 


Table  17.  Variogram  model  parameters  for  elevation  on  a  20  m  grid  for  A.  P.  Hill. 


Variables 

Model  type 

Nugget 

variance 

Sill 

Cl 

Sill 

C2 

Range 
aj  (m) 

Range 

02  (m) 

Raw  elevation 

Pentaspherical 

data 

0 

162.1 

578.0 

Residuals  from 

Double 

linear  trend 

spherical 

0 

19.7 

8.9 

235.0 

927.0 

Residuals  from 

Double 

quadratic  trend 

spherical 

0 

18.7 

14.1 

220.0 

906.0 

Residuals  from 

Double 

Cubic  trend 

spherical 

0 

26.4 

9.5 

263.2 

1363.0 

The  other  aspect  of  these  variograms  is  that  there  is  more  variation  along  the  columns 
than  the  rows.  This  was  also  detected  in  the  variograms  of  the  SPOT  data.  Since  the 
directional  variation  results  in  different  sill  heights  we  did  not  correct  for  it  as  this  kind 
of  anisotropy  (zonal)  requires  stratification.  Figure  72  shows  that  the  ^o  directional 
variograms  overlap  near  to  the  origin  which  is  the  important  point  for  kriging.  It  means 
that  the  kriged  estimates  will  not  be  affected  by  the  directional  difference  at  the  longer 

lag  distances. 


Variance 


Figure  67:  Experimental  variograms  of  elevation  for  the  rows  and  columns  of  a  5  m 
grid,  and  the  average  variogram  of  these  for  A.  P.  Hill. 
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Figure  68:  Experimental  variograms  of  elevation  for  the  rows  and  columns  of  a  20  m 
grid,  and  the  average  variogram  of  these  for  A.  P.  Hill. 


DEM  20m,  Linear  trend. 


Lag  distance  /  m 


Figure  69:  Experimental  variograms  of  elevation  for  the  rows  and  columns  of  a  20  m 
grid,  and  the  average  variogram  of  these  computed  from  the  residuals  of  a  linear  trend, 
and  the  fitted  models  (solid  lines). 


DEM  20m.  Quadratic  trend. 


Figure  70:  Experimental  variograms  of  elevation  for  the  rows  and  columns  of  a  20  m 
grid,  and  the  average  variogram  of  these  computed  from  the  residuals  of  a  quadratic 
trend,  £uid  the  fitted  models  (solid  lines). 
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DEM  20m.  Cubic  trend. 
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Figure  71:  Experimental  variograms  of  elevation  for  the  rows  and  columns  of  a  20  m 
grid,  and  the  average  variogram  of  these  computed  from  the  residuals  of  a  cubic  trend, 
and  the  fitted  models  (solid  lines). 


elevation 


Figure  72:  Experimental  variograms  of  elevation  for  the  rows  and  columns  of  a  20  m 
grid  for  A.  P.  Hill;  D1  is  from  the  rows  and  D2  is  from  the  columns. 
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We  decided  to  do  the  factorial  kriging  on  the  raw  data  and  on  the  residuals  from  the 
quadratic  trend.  Figure  73  shows  the  experimental  variogram  and  the  model  fitted  to  the 
average  variogram  of  the  quadratic  residuals.  Figures  74  and  75  show  the  kriged  maps 
of  the  raw  elevation  data  and  that  computed  using  the  residuals  from  the  quadratic 
trend,  respectively.  They  are  very  similar,  but  the  one  for  the  residuals  shows  slightly 
more  detail.  The  trend  function  that  has  been  removed  varies  very  smoothly  and  this  has 
reduced  the  apparent  variation  in  Figure  73.  The  strange  pattern  in  the  southwestern 
comer  is  the  result  of  missing  data  in  this  region. 

Figures  76  and  77  show  the  maps  of  the  long-range  component  of  the  variation  for  the 
raw  data  and  that  from  the  quadratic  residuals.  The  patterns  that  are  evident  are  similar. 
It  is  interesting  to  compare  these  Figures  with  Figure  80,  the  long  range  component  of 
NIR.  The  large  valley  that  extends  from  the  SW  then  E  and  then  changes  direction  to 
the  N  E  corresponds  with  the  yellow  area  in  Figure  80.  The  tributary  valley  that  extends 
N  from  the  major  valley  at  about  Easting  295500  (Figure  77)  is  also  evident  in  Figure 
80  as  the  orange  area.  The  other  valley  system  to  the  NW  which  mns  in  and  E  to  W 
direction  is  also  evident  in  both  Figures  77  and  80.  The  reddish  areas  on  the  map  of  the 
long  range  component  for  NIR  were  mterpreted  as  higher  ground  in  a  previous  report 
and  this  is  confirmed  now  in  the  maps  of  the  DEM. 

Figures  78  and  79  show  the  maps  of  the  short  range  component  of  the  variation  for  the 
raw  data  and  that  from  the  quadratic  residuals.  The  patterns  are  almost  identical.  These 
correspond  to  the  local  variation  in  relief  that  is  evident  on  the  ordnance  map  for  this 
area.  The  intricacy  of  the  valley  systems  and  general  dissection  is  very  clear  in  these 
maps.  Again  it  is  interesting  to  compare  these  with  the  short  range  component  of  NIR, 
Figure  81.  A  similar  degree  of  detail  and  spatial  scale  is  evident  in  all  of  the  short-range 
maps. 

If  trend  is  present  in  data  it  should  be  removed  for  a  geostatistical  analysis.  In  this  case 
it  has  affected  the  form  of  the  variograms,  but  it  appears  to  have  had  little  effect  on  the 
spatial  patterns  observed  after  kriging  both  raw  and  residual  data. 

The  relation  between  elevation  and  NIR  was  mentioned  above.  This  was  examined 
further  by  computing  the  cross  variogram  between  them.  Figure  80  shows  the 
autovariograms  for  elevation  and  NIR  together  with  their  cross  variogram.  The  outer 
dotted  lines  on  the  cross  variogram  show  the  hull  of  perfect  correlation.  It  shows  that 
although  there  is  a  relation  between  the  image  information  and  elevation  it  is  weak. 
Visually  it  is  more  convincing.  TIus  result  will  be  explored  further  in  the  next  phase  of 
work  when  we  shall  analyse  a  smaller  part  of  the  image. 
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DEM  (20m).  Quadratic  trend. 


Figure  73:  Average  experimental  variogram  of  elevation  for  the  20  m  grid  computed 
from  the  residuals  of  a  quadratic  trend,  and  the  fitted  models  (solid  line). 
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Figure  74:  Kriged  map  of  the  raw  elevation  data  on  the  20  m  grid  for  A.  P.  Hill. 
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Figure  75:  Kriged  map  of  the  residuals  from  a  quadratic  trend  of  elevation  on  the  20  m 
grid  A.  P.  Hill. 
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Figure  76:  Kriged  map  of  the  long  range  component  of  the  variation  for  the  raw 
elevation  data  on  the  20  m  grid  for  A.  P .  Hill. 
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Figure  77:  Kriged  map  of  the  long  range  component  of  the  variation  for  the  residuals 
from  a  quadratic  trend  of  elevation  on  the  20  m  grid  A.  P .  Hill. 
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Figure  78:  Kriged  map  of  the  short  range  component  of  the  variation  for  the  raw 
elevation  data  on  the  20  m  grid  for  A.  P.  Hill. 
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Short  range  component  of  20m  DEM  for  A.  P.  Hill 
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Figure  79:  Kriged  map  of  the  short  range  component  of  the  variation  for  the  residuals 
from  a  quadratic  trend  of  elevation  on  the  20  m  grid  A.  P .  Hill. 
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Figure  80;  Kriged  map  of  the  long  range  component  of  the  variation  for  NIR. 
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Figure  81:  Kriged  map  of  the  short  range  component  of  the  variation  for  NIR, 
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Wavelets  and  Kriging  for  Filtering  and  Data  Reconstruction 


M.  A.  OLIVER*,  E.  BOSCrf  and  K.  SLOCUM^ 

^Department  of  Soil  Science,  The  University  of  Reading,  Whiteknights,  Reading 
RG6  6DW,UK, 

^USTopograhic  Engineering  Center,  7701  Telegraph  Road,  Alexandria,  Virginia 
22310-3864,  U.  S.  A. 

Abstract 

Wavelet  analysis  operates  locally  and  can  describe  a  wide  range  of  frequencies 
simultaneously  and  filter  them  by  multi-resolution  analysis.  Knging  analysis  also  filters 
spatial  variation  at  different  resolutions.  We  compare  the  effectiveness  of  wavelets  and 
factorial  kriging  for  exploring  nested  variation  in  a  SPOT  image.  In  addition  both 
wavelets  and  kriging  can  be  used  to  restore  image  data  after  compression.  We  compare 
the  reliability  of  the  restorations  from  the  two  approaches. 

The  near  infrared  (NDR.)  waveband  of  part  of  a  SPOT  image  covering  Fort  A.  P.  Hill 
in  Virginia  was  used  for  these  analyses.  The  region  is  on  the  dissected  Piedmont  area 
of  the  eastern  United  States.  An  area  of  128  by  128  pixels  was  selected  from  the  scene 
for  analysis.  The  experimental  variogram  was  computed  and  modelled  by  a  nested 
spherical  function  with  correlation  structures  of  about  6.5  pixels  and  21  pixels.  The 
variogram  and  factorial  kriging  separated  the  two  main  spatial  features  present.  The 
low-frequency  component  from  the  wavelet  analysis  contained  the  spatial  structure. 
The  long-range  component  became  evident  as  the  resolution  decreased.  The  high- 
frequency  components  removed  only  the  uncorrelated  variation  and  we  could  not 
retrieve  the  short-range  component. 

The  image  was  sampled  so  that  one  in  every  four  pixels  was  retained,  one  in  every  16 
and  one  in  every  64.  Using  the  variogram  model  for  the  full  set  of  data  values  were 
estimated  at  the  former  data  points  by  ordinary  kriging.  The  low-frequency  wavelet 
transform  for  these  resolutions  was  inverted  so  that  the  missing  values  were  restored. 
The  restored  values  from  both  analyses  were  compared  with  the  original  values  and  the 
mean  squared  differences  (MSD)  computed.  For  all  resolutions  the  MSD  was  smaller 
for  the  wavelet  reconstruction.  However,  the  MSD  proved  somewhat  misleading  when 
frequency  distributions  of  the  errors  were  compared.  They  suggested  that  wavelets  are 
more  able  to  deal  with  the  local  fluctuations  present  in  the  image  and  with  local  non- 
stationarity  than  kriging,  but  that  for  the  majority  of  points  the  kriged  estimates  have  a 
smaller  error. 

The  paper  will  be  illustrated  with  maps  of  the  results,  and  we  shall  suggest 
improvements  for  restoring  images  by  kriging. 


non 


Appendix  III 

♦C  ****  PROGRAM  TO  COMPUTE  VARIOGRAMS  FOR  SQUARES  OF  VARIOUS  SIZES 
C 

C  ****  R  WEBSTER  ROTHAMSTED 
C  Latest  version  22  July  1999 

C 

C  This  program  was  written  as  part  of  US  project 
C  and  may  be  handed  over  to  TEC. 

C 

C  Program  reads  data  on  a  grid  with  X  and  Y  coordinates 

C  and  converts  them  to  an  array  for  the  selected  variate 

C  with  implied  coordinates. 

C 

C  It  tiles  the  grid  into  non -overlapping  squares  of 

C  given  side.  Any  points  to  the  bottom  or  right 

C  of  the  grid  left  over  play  no  role. 

C 

DIMENSION  ZK(190, 189) ,  ZA(30),  GRID(30,30) 

C  ****  ZK(  ,  )  will  hold  grid  of  data. 
character*72  TITLE  (2) 

character*72  INFILE,  0P12,  INll,  FDAT 
data  MAXROW,  MAXCOL/190 , 189/ 
data  IN, INDAT,LP/10, 11, 12/ 

PRINT  *  ,  'WHAT  IS  THE  NAME  OF  THE  STEERING  FILE  ?' 

READ  (5,  '(A)')  INFILE 
OPEN  (INDAT,FILE=INFILE,STATUS='OLD' ) 
print  *  ,  'WHAT  IS  THE  NAME  OF  THE  DATA  FILE  ?  ' 
read  (5, ' (a) ' )  INll 
open  (IN, fiie=INll, status= 'OLD' ) 

PRINT  *, 'WHAT  DO  YOU  WANT  TO  CALL  THE  MAIN  OUTPUT  FILE  ?' 
READ  (5, '(A)')  0P12 
OPEN  (LP,FILE=OP12,STATUS='NEW' ) 

PRINT  *, 'WHAT  DO  YOU  WANT  TO  CALL  THE  SECOND  RESULTS  FILE  ?' 
READ  (5, '(A)')  OP12 
OPEN  (LF8 , FILE=OP12 , STATUS= ' NEW ' ) 

READ  (INDAT,10)  TITLE 
WRITE  (LP,10)  TITLE 
C  WRITE  (LF8,12)  TITLE 

10  FORMAT  (A) 

NVAR  =  int (CYNPUT { INDAT) +0.1) 

NSEL  =  int (CYNPUT (INDAT) +0.1) 

MSIDE  =  int (CYNPUT (INDAT) +0.1) 

MAXLAG=  int (CYNPUT (INDAT) +0.1) 

ZMIS  =  CYNPUT (INDAT) 

ILOG  =  int (CYNPUT (INDAT) +0.1) 
if  (ILOG.eq.l)  SHIFT=CYNPUT (INDAT) 

C  ****  nvAR  is  number  of  variates  in  file. 

C  NSEL  is  the  one  selected  for  analysis. 

C  MSIDE  is  the  side  of  the  square  within  which 

C  averages  are  computed. 

C  MAXLAG  is  the  maximum  lag  distance  of  variograms 

C  ZMIS  is  the  value  used  for  missing  or  blank. 

C  ILOG  =  1  to  transform  to  log  to  base  10. 

C  SHIFT  is  a  value  to  be  added  to  data  to  shift  the  origin 

C  before  taking  logarithms. 

C  ****  Set  data  grid  to  blank 
C  if  (ILOG.eq.l)  ZMIS=logl0 (ZMIS) 

do  20  1=1, MAXROW 
do  20  J=l,MAXCOL 

ZK(I, J) =ZMIS-10000.0 
20  continue 
C 

read  ( INDAT ,10)  FDAT 
C  ****  Read  the  data. 

35  NC=0 


NROW=0 
NCOL=0 
36  NC=NC+1 

read  (IN, FDAT, end=45)  ICOL,  IRON,  (ZA{J) ,  J=1,NVAR) 
if  (ICOL.gt.MAXCOL)  then 
write  (LP,38)  ICOL 
stop 
endif 

if  (IROW.gt .MAXROW)  then 
write  (LP,39)  IROW 
stop 
endif 

38  format  (/lOx, 'ICOL  exceeds  array  bound'//) 

39  format  (/lOx, 'IROW  exceeds  array  bound'//) 
if  (NROW.lt. IROW)  NROW=IROW 
if  (NCOL.lt .ICOL)  NCOL=ICOL 
ZL=ZA(NSEL) 
if  (ILOG.eq.l)  then 

if  (ZL.gt.0.01)  then 
ZL=loglO (ZL+SHIFT) 
else 

ZL=ZMIS-100000 .0 
endif 
endif 

ZK (IROW, ICOL) =ZL 
goto  36 
45  continue 
NC=NC-1 

write  (LP,47)  NC,  NROW,  NCOL 
47  format (//  lOx,  'Number  of  data 

1  lOx,  'Number  of  rows 

2  lOx,  'Number  of  columns 

if  (ILOG.eq.l)  write  (LP,  51)  SHIFT 

51  format  (/lOx, 'DATA  TRANSFORMED  TO  LOG  TO  BASE  10'/ 

1  lOx, 'SHIFT  ',F10.3/) 

ZMAX=-99999999 
ZMIN=999999999 
NN=NC 
ZBAR=0 . 0 
SSQ=0.0 
COUNT=0 . 0 
do  54  1=1, NROW 
do  53  J=l,NCOL 
ZZ=ZK(I,J) 

if  (ZZ.le.ZMIS)  goto  53 

if  (ZMAX.lt.ZZ)  ZMAX=ZZ 
if  (ZMIN.gt.ZZ)  ZMIN=ZZ 
DIF=ZZ-ZBAR 
COUNT=COUNT+l . 0 
ZBAR=ZBAR+DIF/COUNT 
SSQ=SSQ+ (1 . 0-1 . O/COUNT) *DIF*DIF 

53  continue 

54  CONTINUE 
A3=0.0 

do  57  1=1, NROW 
do  56  J=l,NCOL 
ZZ=ZK(I, J) 

if  (ZZ.le.ZMIS)  goto  56 
A3=A3+ (ZZ-ZBAR) **3 

56  continue 

57  CONTINUE 
A2=SSQ/C0UNT 

A3= (A3 /COUNT) / (A2*sqrt (A2) ) 

VAR=SSQ/ (COUNT-1. 0) 


'  ,ilO/ 

'  ,ilO/ 

'  ,ilO/) 


c 

c 


c 


58 


1 

1 

1 

2 

3 
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STD=sqrt (VAR) 

write  (LP,  58)  COUNT,  ZMIN,  ZMAX,  ZBAR,  VAR,  STD,  A3 


format  (/AlOX,  '  Count 

lOx,  '  Minimum  ',fl0.4/ 
lOx,  '  Maximum  ',fl0.4/ 
lOx,  '  Mean  ',fl0.4/ 
lOx,  '  Variance  ',fl2.6/ 
lOx,  '  Standard  deviation  ',fl0.4/ 
lOx,  '  Skewness  ',fl0.4/) 


****  Compute  starting  in  top  left  corner  of  grid. 
NTILER=int (NROW/MSIDE) 

NTILEC=int (NCOL/MSIDE) 


do  300  IR=1,NTILER 
IRS=(IR-1)*MSIDE+1 
IRE=IR*MSIDE 
do  300  IC=1,NTILEC 
ICS=(IC-1)*MSIDE+1 
ICE=IC*MSIDE 


11  =  0 

do  210  I=IRS,IRE 
JJ=0 
11=11+1 

do  210  J=ICS,ICE 
JJ=JJ+1 

GRIDdI,  JJ)=ZK{I,  J) 

210  continue 

Q  ****  Data  are  now  transferred  into  array  GRID(  ,  )  covering 
C  a  small  square  of  side  MSIDE. 

C  Initialize  accumulators.  .  . 

do  220  I=1,MAXLAG 
WLAG(I)=0.0 
GAM (I) =0.0 
WT(I)=0.0 
SUM=0 . 0 
SSQ=0.0 
COUNT=0.0 
220  continue 

do  225  1=1, MSIDE 

do  225  J=l, MSIDE 
ZZ=GRID(I,J) 

if  (ZZ.lt.ZMIS)  goto  225 
COUNT=COUNT+l . 0 
DIF=ZZ-SUM 
SUM=SUM+DIF/COUNT 
SSQ=SSQ+ (1 . 0-1 . O/COUNT) *DIF*DIF 
225  continue 

SSQ=SSQ/ (COUNT- 1 .0) 

SDV=sqrt (SSQ) 

write  (LP,  227)  IRS,  ICS,  SUM,  SSQ,  SDV 
227  format  (//5X,  'CORDINATES  ',  16,16/ 

1  5X,  'MEAN  ',  F12.5/ 

2  5X,  'VARIANCE  ',  F12.5/ 

3  5X,  'ST.  DEVIATION  ',  F12.5/) 

write  (LP,230) 

230  format  (/2x,  'LAG  ANGLE  SEMIVARIANCE  COUNT'/) 
do  255  1=1, MSIDE 
do  255  J=l, MSIDE 
Z1=GRID(I,J) 

if  (Zl.lt.ZMIS)  goto  255 
do  245  K=l, MSIDE 
do  245  L=l, MSIDE 
Z2=GRID(K,L) 


if  (Z2.1t.ZMIS)  goto  245 
Xl=float (I) 

X2=float (K) 

Yl=float (J) 

Y2=float (L) 

D=sqrt ( (X1-X2) **2+(Yl-Y2)**2) 

LAG=int (D) +1 
WLAG  (LAG)  =WIiAG  (LAG)  +D 
GAM ( LAG) =LAG ( LAG) + ( Z1 - Z2 ) * * 2 
WT (LAG) =WT (LAG) +1 . 0 
245  continue 

255  continue 
ANGLE=0 . 0 
do  270  I=1,MAXLAG 

WLAG(I)=WLAG(I)/WT(I) 

GAM (I) =0.5*GAM(I)/WT(I) 

write  (LP,275)  WLAG(I) , ANGLE, GAM(I) ,WT(I) 

270  continue 

275  format  (2x, f 7 . 2 , f 6 . 2 , f 12 . 5, f 10 . 1) 

300  continue 
stop 
end 
C 

FUNCTION  CYNPUT(IN) 

C  ****  READS  A  REAL  NUMBER  FROM  AN  80 -BYTE  RECORD  IN  FREE  FORMAT 
DIMENSION  K(80) ,NUM(10) 

DATA  NUM/lHO , IHl , 1H2 , 1H3 , 1H4 , 1H5 , 1H6 , 1H7 , 1H8 , 1H9/ 

DATA  INOLD,N, IFL,NPLUS,MINUS,NDOT/0, 81, 0, 1H+, IH- , IH. / 
CYNPUT=-0.0 

IFdNOLD.EQ.IN.AND.  N.LE.80)  GOTO  20 
5  IF(IFL.NE.O)  RETURN  .  • 

INOLD=IN 

READ(IN,10)  (K(I) , 1=1, 80) 

10  FORMAT ( 8 OAl) 

15  N=1 

20  IF(N.GT.80)  GOTO  35 
DO  30  I=N,80 
II=K(I) 

DO  25  J=l,10 

IF(II.EQ.NUM(J) )  GOTO  40 
25  CONTINUE 

IF (II. EQ. MINUS)  GOTO  40 
IF(II.EQ.NDOT)  GOTO  40 
IF(II .EQ.NPLUS)  GOTO  40 
30  CONTINUE 
35  GOTO  5 
40  SIGN=1.0 

IF (II .EQ. MINUS)  SIGN=-1.0 
IF  (II. EQ. MINUS  .OR.  II. EQ.NPLUS)  1  =  1  +  1 
IF(I.GT.80)  GOTO  60 
DO  55  N=I,80 
NN=K(N) 

IF(NN.EQ.NDOT)GOTO  70 
DO  45  J=l,10 
KK=J-1 

IF(NN.EQ.NUM(J) )  GOTO  50 
45  CONTINUE 
GOTO  65 

50  CYNPUT=10.0*CYNPUT+KK 
55  CONTINUE 
60  N=82 

65  CYNPUT=SIGN*CYNPUT 
RETURN 
70  I=N+1 


oononnoonno 


TENS=1.0 

IF(I.GT.80)  GOTO  90 
DO  85  N=I,80 
NN=K(N) 

DO  75  J=l,10 
KK=J-1 

IF(NN.EQ.NUM(J) )  GOTO  80 
75  CONTINUE 
GOTO  65 

80  TENS =TENS* 0.1 

CYNPUT=CYNPUT+TENS  *  KK 
85  CONTINUE 
90  N=82 
GOTO  65 
END 

****  PROGRAM  TO  COMPUTE  MOVING  VARIANCES  FOR  SQUARES 
OF  VARIOUS  SIZES 
****  R  WEBSTER  ROTHAMSTED 

Latest  version  22  July  1999 


C 

C 

C 

C 

C 

C 

C 

C 
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This  program  was  written  as  part  of  US  project 
and  may  be  handed  over  to  TEC. 

Program  reads  data  on  a  grid  with  X  and  Y  coordinates 
and  converts  them  to  an  array  for  the  selected  variate 
with  implied  coordinates. 

DIMENSION  ZK(190,189),  ZA(30),  VM(190,189),  AM(190,189) 

****  zK(  ,  )  will  hold  grid  of  data. 
character*72  TITLE (2) 

character* 72  INFILE,  OP12,  INll,  FDAT  .  . 

data  MAXROW,  MAXCOL/190 , 189/ 
data  IN,INDAT,LP/10,11,12/ 

PRINT  *  ,  'WHAT  IS  THE  NAME  OF  THE  STEERING  FILE  ?' 

READ  (5, '(A)')  INFILE 

OPEN  ( INDAT , FILE=INFILE , STATUS= ' OLD ' ) 

print  *  ,  'WHAT  IS  THE  NAME  OF  THE  DATA  FILE  ?  ' 

read  (5, ' (a) ')  INll 

open  (IN,file=INll,status='OLD') 

PRINT  *, 'WHAT  DO  YOU  WANT  TO  CALL  THE  MAIN  OUTPUT  FILE  ? 

READ  (5, ' (A) ' )  OP12 

OPEN  (LP,FILE=OP12,STATUS='NEW' )  r, , 

PRINT  *,  'WHAT  DO  YOU  WANT  TO  CALL  THE  SECOND  RESULTS  FILE 

READ  (5, '(A)')  OP12 

OPEN  (LF8 , FILE=OP12 , STATUS= ' NEW ' ) 

READ  { INDAT ,10)  TITLE 
WRITE  (LP,10)  TITLE 
WRITE  (LF8,12)  TITLE 
10  FORMAT  (A) 

NVAR  =  int(CYNPUT( INDAT) +0.1) 

NSEL  =  int(CYNPUT( INDAT) +0.1) 

MSIDE  =  int(CYNPUT( INDAT) +0.1) 

ZMIS  =  CYNPUT (INDAT) 

ILOG  =  int (CYNPUT (INDAT) +0.1) 
if  (ILOG.eq.l)  SHI FT=CYNPUT (INDAT) 

****  NVAR  is  number  of  variates  in  file. 

NSEL  is  the  one  selected  for  analysis. 

MSIDE  is  the  side  of  the  square  within  which 
averages  are  computed. 

ZMIS  is  the  value  used  for  missing  or  blank. 

ILOG  =  1  to  transform  to  log  to  base  10.  _ 

SHIFT  is  a  value  to  be  added  to  data  to  shift  the  orig 
before  taking  logarithms. 

****  Set  data  grid  to  blank 


